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FOREWORD 


NASTRAN®  (NASA  STRUCTURAL  ANALYSIS)  is  a  large,  comprehensive, 
nonproprietary,  general  purpose  finite  element  computer  code  for  structural 
analysis  which  was  developed  under  NASA  sponsorship  and  became  available  to 
the  public  in  late  1970.  It  can  be  obtained  through  COSMIC®  (Computer 
Software  Management  and  Information  Center),  Athens,  Georgia,  and  is  widely 
used  by  NASA,  other  government  agencies,  and  industry.^ 

NASA  currently  provides  continuing  maintenance  of  NASTRAN  through  COSMIC. 
Because  of  the  widespread  interest  in  NASTRAN,  and  finite  element  methods  in 
general, \the  Eighteenth  NASTRAN  Users'  Colloquium  was  organized  and  held  at 
The  Red.  Li'on  Portland  Center,  Portland,  Oregon  on  April  23-27,  1990.  (Papers 
from  previob-s^  colloquia  held  in  1971,  1972,  1973,  1975,  1976,  1977,  1978, 
1979,  1980,  19820.983,  1984,  1985,  1986,  1987,  1988  and  1989  are  published  in 
NASA  Technical  Memorandums  X-2378,  X-2637,  X-2893,  X-3278,  X-3428,  and  NASA 
Conference  Publication  2018,  2062,  2131,  2151,  2249,  2284,  2328,  2373,  2419, 
2481,  2505  and  3029.)  Mhe  Eighteenth  Colloquium  .provides  some  comprehensive 
general  papers  on  the  application  of  finite  element  methods  in  engineering, 
comparisons  with  other  approaches,  unique  applications,  pre-  and 
post-processing  or  auxiliary  programs,  and  new  methods  of  analysis  with 


NASTRAN.  — — ^  Q>1‘  ^ 

Individuals  actively  engaged  in  the  use  of  'inite  elements  or  NASTRAN 
were  invited  to  prepare  papers  for  presentation  at  the  Colloquium.  These 
papers  are  included  in  this  volume.  No  editorial  review  was  provided  by  NASA 
or  COSMIC;  however,  detailed  instructions  were  provided  each  author  to  achieve 
reasonably  consistent  paper  format  and  content.  The  opinions  and  data 
presented  are  the  sole  responsibility  of  the  authors  and  their  respective 
organizations. 


NASTRAN®  and  COSMIC®  are  registered  trademarks  of  the  National  Aeronautics  and 
Space  Administration. 
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A  NASTRAN  TRAINER  FOR  DYNAMICS 

H.R.  Grooms,  P.J.  Hinz,  and  G.L.  Commerford 
Rockwell  International 


OVERVIEW 

As  the  use  of  the  finite  element  method  proliferates,  the  need  for  training 
becomes  more  and  more  pronounced.  An  automated  tool  to  familiarize  engineers  with 
static  solutions  has  been  developed  and  used.  This  tool  (Figure  1)  is  part  of  ah 
overall  structural  analysis/expert  training  system  (ref.  1).  Experiences  with  this 
tool  and  comments  from  users  (ref.  2).  have  underlined  the  need  for  a  dynamic  version 
of  the  trainer. 

This  paper  presents  an  automated  training  tool  that  engineers  can  use  to  master 
the  application  of  NASTRAN  to  dynamic  problems.  The  paper  consists  of  the  following 
sections: 

•  Overview 

•  Background 

•  Existing  Programs 

•  Scope,  Purpose,  and  System  Organization 

•  Example  Problems 

•  Conclusions 

•  References 

Example  problems  have  been  selected  to  make  classical  solutions  available  for 
comparison.  These  comparisons  can  be  used  to  evaluate  the  solution. 


BACKGROUND 

The  solution  of  dynamic  problems  involves  some  complications  that  do  not  exist 
with  static  problems: 

•  How  many  degrees  of  freedom  should  be  retained  for  the  eigenvalue  solution? 

•  Which  discrete  mass  items  are  so  large  or  important  that  they  should  be 
retained  for  eigenvalue  solution? 

•  How  many  frequencies  and  mode  shapes  are  needed  and  to  what  accuracy? 

An  engineer  may  think  that  most  of  the  mass  associated  with  a  structure  can  be 
traced  to  the  structural  members  themselves;  this  is  not  necessarily  true.  With  many 
aircraft  and  spacecraft,  the  nonstructural  masses  (e.g.,  hydraulic,  lines,  fuel  tanks, 
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environmental  control  equipment,  etc.)  have  a  pronounced  influence  on  the  overall 
mass  distribution  and  may  have  the  greatest  dynamic  effect. 

The  example  problems  have  distributed  masses  and  lumped  masses  that  the  user 
must  consider  in  the  solution  approach.  These  examples  help  the  user  develop  judgment 
when  deciding  on  the  number  and  the  particular  degrees  of  freedom  to  be  retained,  and 
on  how  to  discretize  the  distributed  mass. 


EXISTING  PROGRAMS  • 

! 

Various  researchers  have  developed  computer  programs  for  structural  analysis  and 
design  applications.  Ginsburg  (ref.  3)  'addressed  computer  literacy,  while  Woodward 
and  Morris  discussed  improved  productivity  through  interactive  processing  (ref.  4). 
Wilson  and  Holt  (ref.  5)  developed  a  system  for  computer-assisted  learning  in 
structural  engineering.  Sadd  and  Rolph  (ref.  6)  described  the  various  ways  in  which 
design  engineers  could  be  trained  to  use  the  finite  element  method.  Self-adapting 
menus  for  CAD  software  are  covered  by  Ginsburg  (ref.  7). 

Bykat  (ref.  8)  is  developing  a  system  that  will  have  features  for  training, 
analysis  control,  and  interrogation. 


SCOPE,  PURPOSE,  AND  SYSTEM  ORGANIZATION 

The  NASTRAN  trainer  was  designed  to  be  a  stand-alone  tool.  The  trainer  is  user 
friendly — a  knowledge  of  job  control  language  or  the  operating  system  is  not 
required.  A  user  can  sit  down  at  a  terminal  and,  in  very  little  time,  start  solving 
an  example  problem.  The  trainer  is  organized  so  that  a  user  must  complete  the  static 
problems  (ref.  2)  before  the  dynamics  problems  can  be  accessed.  This  organization 
prevents  a  user  who  has  no  familiarity  with  the  finite  element  method  from  starting 
with  the  dynamics  section. 

The  trainer  is  organized  into  three  main  modules:  (1)  overview,  (2)  user's 
guide,  and  (3)  problem  set.  Figure  2  shows  some  details  of  each  module.  The  user 
accesses  these  modules  by  using  the  primary  menu.  More  details  of  the  NASTRAN 
environment  sections  are  given  in  Figure  3. 


EXAMPLE  PROBLEMS 

The  example  problems,  shown  in  Figures  4  through  11  and  summarized  in  Table  1, 
become  progressively  more  difficult  to  solve.  The  first  problem  is  a  simply  supported 
beam  with  a  single  lumped  mass  at  the  center. 

There  are  various  courses  and  classes  to  instruct  engineers  in  solving  dynamics 
problems.  These  courses  usually  emphasize  the  theory.  A  vital  part  of  solving  any 
large  dynamics  problems  is  deciding  how  many  and  which  degrees  of  freedom  should  be 
retained  for  the  eigenvalue  solution.  This  is  usually  a  matter  of  judgment,  and  it 
takes  solving  many  problems  to  develop  this  judgment. 
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Example  2  was  solved  using  three  different  approaches.  The  user  was  trying  to 
answer  some  fundamental  questions  that  must  be  addressed  every  time  a  dynamics 
problem  is  solved  using  the  finite  element  method: 

•  Is  the  model  fine  enough? 

•  Have  the  distributed  masses  been  lumped  into  enough  locations? 

•  Have  enough  degrees  of  freedom  been  retained  in  the  eigensolution? 

Figure  12  summarizes,  the  different  approaches.  Table  2  compares  the  computed 
three  lowest  natural  frequencies  with  the  exact  results. 


CONCLUSIONS 

An  automated  training  tool  that  helps  engineers  become  familiar  with  using 
NASTRAN  to  solve  dynamic  problems  has  been  presented.  The  tool  allows  the  user  to 
proceed  at  his  own  pace  by  using  a  set  of  eight  example  problems.  The  examples  were 
selected  so  that  classical  solutions  are  available  and  displayed,  enabling  the  user 
to  make  comparisons. 
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Table  1.  Example  Problems 


Example 

Description 

Significant  Features 

1 

Beam  simply  supported  on  both  ends  with  lumped  mass  in  middle 

Motion  in  one  plane  only,  lumped  mass  only 

2 

Beam  simply  supported  on  both  ends  with  uniformly  distributed 
mass 

Motion  in  one  plane  only,  distributed  mass 

3 

Beam  fixed  on  one  end  with  a  lumped  mass  at  the  free  end 

Motion  in  any.direction,  lumped  mass  only 

4 

Beam  fixed  on  one  end  with  a  uniformly  distributed  mass 

Motion  in  any  direction  with  uniformly  distributed  mass* 

5 

Rectangular  plate  clamped  on  one  edge,  all  other  edges  free  with 
a  uniformly  distributed  mass 

Plate  bending  with  distributed  mass 

6 

Rectangular  plate,  free-free  with  uniformly  distributed  mass 

Free-free  (implies  six  modes  with  zero  frequency) 

7 

Two  beams  connected  by  springs,  each  with  distributed  and 
lumped  mass 

Multibody  problem,  free-free 

8 

Problem  7  with  a  forcing  function  added 

Forcing  function 

Table  2.  Comparison  of  Natural  Frequencies  for  Example  2 


”  ^^Approach 
Mode 

First  Approach 

Second 

Third 

1 

9.870 

9.867 

9.869 

9.872 

2 

39.48 

39.19 

39.47 

39.74 

3 

88.83 

83.21 

88.66 

93.62 

i 


*  ' 
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Figure  2.  Organization  of  NASTRAN  Trainer 
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Figure  3. 


Organization  of  Program 


Figure  4.  Simply  Supported  Beam  with  Concentrated  Mass  (Example  1) 


w=  LB/FT 


CROSS  SECTION 


Figure  5.  Simply  Supported  Beam  with  Uniformly  Distributed  Mass  (Example  2) 
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CROSS  SECTION 


Figure  6.  Cantilever  Beam  with  Concentrated  Mass  (Example  3) 


CROSS  SECTION 


Figure  7.  Cantilever  Beam  with  Uniformly  Distributed  Mass  (Example  4) 


A)  E,  1 1  , 1  2 


Li  Li  Li 


Figure  10.  Two  Beams  Connected  by  Two  Springs  (Example  7) 
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Figure  11.  Two  Beams  Connected  by  Two  Springs  Driven 
by  a  Forcing  Function  (Example  8) 
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DOFs  IN 

APPROACH 

NO.  OF 

NO.  OF 

EIGEN- 

ELEMENTS 

MASSES 

SOLUTION 

wL/8 

wL/4  wL74 

wL/4 

wU8 

1 

IHzl 

■ 

4 

5 

3 

( 

U4  L/2 

3L/4 

L 

L/8  3L78 

5L/8 

7L78 

wL/8  EACH 

wL/16 

1 

I  wL/16 

2 

■ 

8 

9 

7 

■ 

3 

SAME  AS  APPROACH  2,  EXCEPT  SAME  DOFs  SAVED 

8 

9 

3 

FOR  EIGENSOLUTiON  AS  APPROACH  1 

Figure  12.  Three  Approaches  to  Example  2 
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PCI  -  A  PATRAN-NASTRAN  MODEL  TRANSLATOR 


T.J.  SHEERER 

CHRYSLER  TECHNOLOGIES  AIRBORNE  SYSTEMS,  INC. 

WACO,  TX 


INTRODUCTION 


The  existence  of  several  derivative  versions  of  NASTRAN,  which 
differ  significantly  in  element  definitions  and  result  formulation, 
has  caused  some  difficulties  in  the  interface  between  NASTRAN  and 
pre-  and  post-  processors  such  as  PATRAN  or  SUPERTAB.  In 
particular,  the  PATRAN-COSMIC/NASTRAN  interface  provided  by  PDA 
Engineering  has  not  been  updated  at  the  same  rate  as  the  equivalent 
interface  with  MSC/NASTRAN,  and  has  significantly  less 
capabilities.  Many  model  entities  supported  by  both  PATRAN  and 
COSMIC/NASTRAN  are  not  supported  by  the  translator.  The  well- 
documented  PATRAN  neutral  file,  which  is  now  supported  by  several 
other  vendors,  has  provided  a  means  for  the  user  to  create  his  own 
interface  program  for  model  translation,  while  it  has  also  been 
possible  to  pass  results  form  NASTRAN  to  PATRAN  with  a  user-written 
program  using  0UTPUT2  statements  and  format  information  from  the 
Patran  Users  1  Manual .  In  recent  years  PDA  engineering  has  provided 
as  part  of  PATRAN  a  library  of  subroutines  known  as  gateway 
utilities,  which  extract  data  directly  from  the  PATRAN  database 
file  and  which  can  be  called  from  a  FORTRAN  program.  As  this 
eliminates  the  task  of  reading  the  neutral  file  the  work  of 
creating  a  translator  can  be  produced  by  the  user  with  little 
effort.  This  has  been  done  with  the  object  of  producing  a  PATRAN- 
COSMIC  NASTRAN  model  translator  comparable  in  scope  to  the  PATRAN- 
MSC/NASTRAN  translator,  but  also  allowing  a  greater  degree  of  user 
control  than  is  found  therein.  The  different  parts  of  the  program 
were  developed  in  several  locations,  as  the  counterpart  to  a 
results  translator  developed  for  Texas  Instruments  by  Texas  A  and 
M  University,  which  also  emphasizes  flexibility.  Both  these 
programs  are  public  domain  programs  under  the  terms  of  the 
development  agreement  between  Texas  Instruments  and  Texas  A.M.U., 
and  enhancements  developed  at  Chrysler  have  also  been  passed  to 
T.A.M.U. 
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PROGRAM  STRUCTURE 


PCI  supports  a  range  of  elements  comparable  with  PATNAS  and 
significantly  greater  that  PATCOS*.  The  structure  of  the  program 
is  such  that,  'in  effect,  it  supports  elements  not  currently  extant. 
In  the  PATCOS  program  element  types  are  hard-wired,  and  if  a 
different  NASTRAN  card  is  required  a  text  editor  must  be  used  to 
change  the  bulk  data  file.  If,  for  example,  the  QUAD4  element  is 
required  in  a  model,  all  elements  of  this  type  must  be  edited  from 
type  QUAD2  to  QUAD4,  since  PATCOS  does  not  currently  support  the 
QUAD 4 .  Similar  difficulties  arise  with  other  elements.  The  PATNAS 
translator  does  support  a  wide  range  of  elements  identified  by 
configuration  codes,  but  the  equivalence  between  elements  and 
configuration  code  is  controlled  from  within  PATRAN,  and  the  user 
can  not  translate  to  an  element  developed  in-house,  or  a  newly 
introduced  element,  or  simply  an  element  not  supported  by  the 
translator.  PCI  uses  a  different  approach  by  having  a  user  text 
file  of  twenty  NASTRAN  element  names  for  each  element 
configuration.  This  file  may  be  edited  by  the  user  to  assign 
whatever  correspondence  he  desires  between  configuration  code  and 
NASTRAN  element  type.  Since  most  element  cards  in  NASTRAN  follow 
the  same  pattern  of  (NAME, PID/MID, NODES)  it  is  not  generally 
necessary  to  write  a  new  subroutine  when  adding  an  element. 
Frequently  the  only  parameter  changed  is  field  one  of  the  NASTRAN 
data  card,  and  this  is  controlled  from  a  text  file.  Because 
NASTRAN  and  PATRAN  use  a  different  numbering  sequence  for  higher 
elements  it  is  convenient  to  use  several  subroutines  to  write 
NASTRAN  elements,  but  frequently  only  one  routine  is  required  for 
each  shape/node  combination.  All  3-node  shells,  for  example,  are 
written  by  the  same  routine.  The  PATRAN  configuration  field  is 
used  to  select  an  element  name  from  a  range  of  twenty  for  each 
configuration  of  shape/nodes.  The  user-generated  text  file 
ETYPES.DAT  contains  these  names  for  configurations  1  through  20. 
Any  elements  having  the  default  configuration  of  zero  in  the  PATRAN 
database  are  assigned  the  value  1  in  the  translator. 

Table  (1)  shows  the  basic  PATRAN  element  types  supported  by 
PCI  and  the  NASTRAN  elements  obtainable  form  them.  For  comparison 
the  elements  supported  by  PATCOS  and  PATNAS  are  also  shown,  note 
that,  in  general,  failure  by  PCI  to  support  elements  listed  for 
PATNAS  is  because  they  are  not  supported  in  COSMIC/NASTRAN.  In 
practice,  alteration  of  the  text  file  of  element  configurations 
will  allow  support  of  these  elements  if  and  when  they  become 
available 
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TABLE  1 :  GEOMETRIC  ELEMENT  TRANSLATION 

CAPABILITIES 

OF  PATCOS,  PATNAS  AND  PCI 

SHAPE/NODES 

PATCOS 

PATNAS 

PCI 

BAR/ 2 

CBAR 

CBAR 

CBAR 

CBEAM 

CROD 

CROD 

TRIA/3 

CTRIA2 

CTRIAl 

CTRIAl 

CTRIA2 

CTRIA3 

CTRIA2 

CTRBSC 

CTRIARG 

CTRIARG 

CTRMEM 

CTRPLT 

QUAD/ 4 

CQUAD2 

CQUAD1 

CQUAD1 

CQUAD2 

CQUAD2 

CQUAD4 

CQUAD4 

CQDMEM1 

CQDMEM1 

CQDMEM2 

CQDMEM2 

CSHEAR 

CSHEAR 

CTWIST 

CQDPLT 

CQDMEM 

CTRAPARG 

CTRAPAX 

TRIA/6 

CTRIM6 

CTRIA6 

CTRIM6 

CTRIAX6 

CTRSHL 

CTRPLT1 

QUAD/ 8 

CIS2D8 

CQUAD8 

CIS2D8 

HEX/ 8 

CIHEX1 

CHEXA 

CIHEX1 

CHEX8 

CHEXA1 

CHEXA2 

QUAD/ 20 

CIHEX2 

CHEXA 

CIHEX1 

CHEX8 

CHEXA1 

CHEXA2 

QUAD/ 3 2 


CIHEX3 
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PCI  also  supports  the  SCALAR  and  damping  elements  CELAS2  and 
CDAMP2  by  generating  a  scalar  element  for  each  degree  of  freedom 
specified  in  a  PATRAN  SPRING  element.  The  concentrated  mass 
element  C0NM2  is  obtainable  from  the  PATRAN  MASS  directive.  MPCs 
and  rigid  elements  are  supported.  Table  (2)  summarizes  the  support 
for  these  elements.  Node  translation  with  embedded  SPCs  is  fully 
supported . 


TABLE  2: 

SCALAR,  DAMPER  AND  MASS 

ELEMENT  SUPPORT 

PATRAN  DIRECTIVE 

NATRAN 

CARD  WRITTEN  BY  PCI 

MPC 

MPC 

MPC , RROD 

CRIGIDR 

MPC , RBAR 

CRIGD1 

MPC, RBE1 

CRBE1 

MPC, RBE2 

CRBE2 

MPC, RBE3 

CRBE3 

BAR/2/n  OR 

SPRING 

CEL AS 2 

BAR/2/n  OR 

MASS 

C0NM2 

BAR/2/n 

CDAMP2 

LOAD  AND  CONSTRAINT  TRANSLATION: 

PCI  supports  constraint  (SPC),  nodal  force  and  pressure 
translation.  FORCE  and  SPC  cards  are  translated  in  an  element- 
dependent  manner  as  shown  in  Table  ( 3 ) .  Only  normal  element 
pressures  are  supported.  PCI  will  select  the  appropriate  PLOAD 
card  for  the  element  type.  In  many  circumstances  involving  higher 
order  elements  several  PLOAD  cards  must  be  generated  for  a  single 
PATRAN  pressure  load. 
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TABLE  3t, 

PCI  SUPPORT  OF  LO^DS  and  spcs 

NASTRAN  CARD 

USED  IN  ELEMENTS 

■SUPPORT 

PL0AD4- 
P LOAD 2 

PL0AD3 

PLOAD 

CQUAD4 

OTHER  1ST  ORDER  SHELL 
ISOPARAMETRIC  SOLID 

2ND  ORDER  SHELL 

■  ~  ,  .. .  _ , _  ,,  ... 

YES 

YES 

NOT  YET 

YES 

COORDINATE  SYSTEMS: 


The  various  PATRMf  coordinate  .systems  are  translated  to  C0RD2 
cards  in  NASTRAN,  as  in  the  PATNAS  translator.  Additionally, 
PATRAN  nodes,  which  are  stored  in  the  database  in  the  basic 
coordinate  system,  have  their  locations  output  to  NASTRAN  in  the 
local  system  associated  with  the  nodes.  This  is  of  great 
importance  if,  for  example,  constraints  are  to  be  applied  in  a 
local  system. 

The  PATRAN  database  includes „  associated  with  each  coordinate 
system,  a  3x3  matrix  T  such  that 


where  the  suffixes  1,  b,  denote  local  and  basic  coordinate  systems 
respectively.  The  translator  inverts  the  matrix  to  produce  the 
matrix  T~'  whore 

[T 


and  prints  the  local  values  to  the  NASTRAN  data  deck. 


PROPERTY  GENERATION: 

It  is  generally  no  more  difficult  to  enter  property  and 
material  cards  in  NASTRAN  that  in  PATRAN.  For  this  reason  property 
generation  has  not  been  implemented  in  PCI. 
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CONCLUSIONS: 


-  The  .'amount  of  programming  r  paired .'.to  develop  a  PATRAN-NASTRAN 
translator  Was  surprisingly  sigai.1 .  The  .approach,  taken  produced  a 
highly  flexible  translator  comparable  -with  the  PATNA3  translator, 
and  superior  to  the  PATCOS  translator.  The  coding  required  varied 
from  around  ten  lines  for  a  shell  element  to  arotnd;  thirty  for  a 
bar  element,  and  the  time  required  to  add  a  feature  to  the  program 
is  typically  less  than  an  hour.  The  use  of  a  lookup  table  for 
element  names  makes  the  translator  also  applicable  to  ether 
versions  of  NASTRAN.  The  saving  in  time  as;  a  result  of  usitvg  PDA's 
Gateway  utilities  was  considerable. 

During  the  writing  of  the  program  it  became  apparent  that, 
with  a  somewhat  more  complex  structure,  it  would  be  possible  to 
extend  the  element  data  file  to  contain  all  data  required  to  define 
the  translation  from  PATRAN  to  8ASTR&&  by  mapping  of  data  between 
formats.  Similar  data  files  on  property,  material  end  grid  formats 
would  produce  a  completely  universal  translator  from  PA TRAN  to  any 
FEA  program,  or  indeed  any  CAE  system. 
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A  GENERIC  iNTERFACE  BETWEEN  COSMIC/NASTRAN  AND  PATRAN 


Paul  N.  Roschke,  Prakit  Premthamkorn,  and  James  C.  Maxwell 
Texas  A&M  University 

SUMMARY 

Despite  its  powerful  analytical  capabilities,  COSMIC/NASTRAN  lacks  adequate 
post-processing  adroitness.  PATRAN1,  on  the  other  hand  is  widely  accepted  for  its 
graphical  capabilities.  A  nonproprietary,  public  domain  code  mnemonically  titled  CPI  (for 
COSMIC/NASTRAN-PATRAN  Interface)  is  designed  to  manipulate  a  large  number  of 
files  rapidly  and  efficiently  between  the  two  parent  codes.  In  addition  to  PATRAN's  results 
file  preparation,  CPI  also  prepares  PATRAN's  P/PLOT  data  files  for  xy  plotting.  The  user 
is  prompted  for  necessary  information  during  an  interactive  session.  Current 
implementation  supports  NASTRAN's  displacement  approach  including  the  following  rigid 
formats:  (1)  static  analysis,  (2)  normal  modal  analysis,  (3)  direct  transient  response,  and 
(4)  modal  transient  response.  A  wide  variety  of  data  blocks  are  also  supported.  Error 
trapping  is  given  special  consideration.  A  sample  session  with  CPI  illustrates  its  simplicity 
and  ease  of  use. 


INTRODUCTION 

Overview 

The  standard  gateway  that  interfaces  COSMIC/NASTRAN's  analysis  results  to 
PATRAN's  post-processing  makes  use  of  NASTRAN's  FORTRAN-written  results  files. 
These  files  can  be  requested  via  DMAP  ALTER's  OUTPUT2  statement  in  NASTRAN's 
executive  data  deck.  They  contain  data  in  mixed  ASCII  and  binary  code  format.  However, 
they  can  not  be  used  as  direct  input  to  PATRAN.  Similarly,  PATRAN  also  supports 
communications  with  external  codes  via  specially  formatted  results  files.  Format  of  these 
files  is  predetermined  according  to  PATRAN  and  differs  for  each  data  type.  Generally, 
they  can  be  categorized  into  three  groups  according  to  their  formats:  (1)  nodal  results  files, 
(2)  element  results  files,  and  (3)  beam  results  files.  The  number  of  results  files  can  be  as 
many  as  required.  Therefore,  in  order  to  interface  COSMIC/NASTRAN  to  PATRAN  for 
post-processing  purposes,  an  interface  that  is  capable  of  translating  from  NASTRAN's 
analysis  results  to  PATRAN-recognizable  files  is  required. 

Fig.  1  shows  code  and  file  relationships  among  the  analysis  and  post-processing 
modules.  The  flow  begins  with  the  input  of  NASTRAN's  analysis  data  deck.  DMAP 
ALTER  statements  must  be  included  in  the  executive  control  deck  in  order  to  obtain  a 
FORTRAN-written  results  file.  Generally,  this  file  is  assigned  an  extension  of  .PAT.  This 
NASTRAN  results  file  is  then  translated  by  the  COSMIC/NASTRAN-PATRAN  Interface 
(CPI).  Output  from  CPI  is  a  group  of  PATRAN-compatible  results  files  and/or  a 
P/PLOT-compatible  data  file.  These  files  are  ready  for  PATRAN  post-processing. 


1  PATRAN  is  a  trademark  of  PDA  Engineering 
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PATRAN  also  requires  access  to  a  neutral  file  containing  geometrical  properties  of  the 
model.  This  file  can  be  obtained  via  the  COSPAT2  translator,  but  not  by  means  of  CPI. 


Fig.  1  Code  and  File  Relationships 


Program  Description 

COSMIC/NASTRAN-PATRAN  Interface  (CPI)  translates  COSMIC/NASTRAN's 
FORTRAN-written  results  files  to  PATRAN  compatible  results  files.  These  results  files 
can  be  requested  via  the  OUTPUT2  instruction  in  the  executive  control  deck.  CPI 
provides  the  operation  in  two  modes: 

Mode  1  -  Whole  Model  Translation:  CPI  translates  all  data  blocks  contained  in  the 
NASTRAN  results  file  (*.PAT)  and  creates  PATRAN  compatible  results  files 
corresponding  to  data  blocks  found  in  the  input  data  file.  The  user  is  prompted  for  the 
prefix  name  of  an  output  file.  Different  prefixes  allow  the  user  to  distinguish  between 
groups  of  output  files  when  many  results  files  are  translated  in  a  single  execution.  Up  to  6 
characters  per  prefix  is  acceptable.  Created  output  files  and  translated  data  blocks  are 
summarized  on  the  screen  during  CPI's  execution. 

Mode  2  -  XY-PIot  Data  Preparation:  CPI  creates  data  for  the  specified  node  or  element 
from  a  set  of  data  blocks.  CPI  prompts  the  user  for  necessaiy  information,  i.e.  node 
number  and  component  number  for  nodal  data,  and  searches  the  input  file  for  any  data  on 


2  COSPAT  is  a  COSMIC/NASTRAN  to  PATRAN  interface,  developed  by  COSMIC,  University  of  Georgia, 
Athens,  Georgia. 
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this  node.  Each  data  set  is  written  to  PATRAN's  P/PLOT  compatible  data  file.  CPI 
names  this  file  XYPLOT.DAT.  When  this  command  is  selected,  CPI  automatically  reviews 
the  data  blocks  available  in  the  current  input  file.  The  user  is  then  prompted  to  enter  the 
number  corresponding  to  the  data  block  name  shown  above  this  prompt. 

Currently,  CPI  supports  NASTRAN’s  displacement  approach  including,  but  not 
limited  to,  the  following  rigid  formats: 

(1)  Static  Analysis  (Rigid  Format  1) 

(2)  Normal  Modal  Analysis  (Rigid  Format  3) 

(3)  Direct  Transient  Response  (Rigid  Format  9) 

(4)  Modal  Transient  Response  (Rigid  Format  12) 

Table  1  shows  data  blocks  supported  by  CPI.  It  should  be  noted  that  CPI  supports 
any  rigid  format  as  long  as  the  data  blocks  listed  are  encountered.  Rigid  formats  named 
above  are  only  a  guide  to  the  rigid  formats  most  frequently  giving  rise  to  these  data  blocks. 
CPI  recognizes  data  blocks,  not  rigid  formats.  Table  2  shows  elements  currently  supported 
by  CPI. 


Table  1.  Data  Blocks  Supported  by  CPI 


Data  Block 

Content 

OUGVl 

Nodal  Displacements 

HOUGV1 

Nodal  Temperature 

OQG1 

Single  Point  Constraint 

OPG1 

Load  Vectors 

OPHIG 

Eigenvectors 

OUPV1 

Displacement  Vectors 

ONRGY1 

Strain  Energy 

OESl(X) 

Element  Stress 

OEF1 

Element  Forces 

OGPFB1 

Grid  Point  Balance  Forces 

Table  2.  Elements  Supported  by  CPI 


Element  Name 

ID  Number 

Element  Name 

ID  Number 

CROD 

l 

CQDPLT 

15 

CQDMEM2 

63 

CBEAM 

2 

CTRIA2 

17 

CTUBE 

3 

CQUAD2 

18 

CTRIA1 

6 

CQUAD1 

19 

CTRBSC 

7 

CBAR 

34 

CTRPLT 

8 

CTRIARG 

36 

CTRMEM 

9 

CTRAPRG 

37 

CONROD 

10 

CQDMEM1 

62 

CTETRA 

39 

CWEDGE 

40 

CHEXA1 

41 

CHEXA2 

42 

CIHEX1 

65 

CIHEX2 

66 

CIHEX3 

67 
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PROGRAM  USAGE  AND  SAMPLE  SESSION 

COSMIC/NASTRAN-PATRAN  INTERFACE  (CPI)  is  an  interactive  program.  A 
user  can  manipulate  the  program  by  selecting  from  menu  commands  and  answering 
questions.  No  special  commands  are  needed.  CPI  also  provides  extensive  error  trapping  to 
ensure  appropriate  input  and  output. 

The  following  sample  session  (adapted  from  Ref.  5)  illustrates  execution  of  CPI  for 
transient  dynamic  analysis  results  in  which  a  series  of  time  steps  is  involved.  A  simple  three 
node  model  which  is  350  cm.  high,  and  has  a  cross  sectional  area  of  6.45  cm.2  (see  Fig.  2),  is 
used  to  model  a  rocket  trajectory.  Time  versus  loading  of  the  rocket  is  shown  in  Fig.  3. 

node 


Fig.  2  Sample  2  -  Geometric  Properties 


Load  (N) 


After  invocation,  CPI's  opening  header  appears  on  the  screen  and  CPI  prompts  the 
user  for  the  name  of  the  NASTRAN  results  file  to  be  translated.  The  user  then  enters  the 
NASTRAN  results  file  name  with  its  extension.  Alternatively,  the  user  can  leave  the 
translator  and  return  to  the  operating  system  by  entering  "EXIT1  (or  "exit"  or  "E"  or  simply 
"e").  Here,  the  NASTRAN  results  file  name  is  assigned  the  name  "ROCKET.PAT." 
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****************************************************** 

** 

** 

** 

COSMIC/NASTRAN  -  PATRAN  INTERFACE 

** 

** 

** 

** 

**  CPI** 

** 

** 

** 

** 

TEXAS  ASM  UNIVERSITY 

** 

** 

& 

** 

** 

TEXAS  INSTRUMENTS 

** 

** 

iri f 

** 

(JUNE  1989) 

** 

** 

** 

****************************************************** 

ENTER  NASTRAN  RESULTS  FILE  NAME  OR  ''EXIT" 

>ROCKET . PAT 

Once  the  results  file  has  been  named,  CPI  reads  the  results  file  header  and  echoes 
it  to  the  screen  for  verification.  Next  a  menu  and  a  prompt  for  selection  appear  on  the 
screen.  The  response  to  CPI  should  be  one  of  the  following: 

>  1  Select  1  to  translate  the  entire  file  and  create  appropriate  PATRAN-compatible 
results  files. 

>  2  Select  2  to  translate  the  data  and  create  PATRAN's  P/PLOT-compatible  data  files 
which  contain  requested  nodal  or  element  data. 

>  3  Select  3  to  specify  a  new  NASTRAN  results  file. 

>  4  Select  4  to  exit  CPI  and  return  to  the  operating  system. 

If  the  first  option  is  chosen,  translation  proceeds  as  follows: 


DATE  : 

7/27/89 

HEADER: 

NASTRAN  FORT  TAPE  ID  CODE  - 

LABEL  : 

xxxxxxxx 

>»»>> 

ENTER  YOUR  SELECTION  <««« 

1) 

WHOLE  MODEL  TRANSLATOR 

2) 

XY-PLOT  DATA 

3) 

ENTER  NEW  RESULTS  FILE 

4) 

EXIT 

>  i 

ENTER  FILENAME  PREFIX  (UP  TO  6  CHARACTERS)  OR  [RETURN] 
(EXAMPLE  DEFAULT  FORMAT:  S  I  .DIS) 

>  ROCKET 

IS  "ROCKET"  THE  CORRECT  PREFIX?  [Y] 

><RETURN> 

1. 


FILE  CREATED 

FROM  DATA  BLOCK 

ROCKETS01I1.DIS 

OUPV1 

ROCKETS01I2 .DIS 

OUPV1 

ROCKETS01I3.DIS 

OUPV1 

II 

II 

It 

II 

It 

II 

ROCKETSO 1150. DIS 

OUPV1 

ROCKETSO 1151. DIS 

OUPV1 

»>»»  ENTER  YOUR  SELECTION  <««« 


1)  WHOLE  MODEL  TRANSLATOR 

2)  XY-PLOT  DATA 

3)  ENTER  NEW  RESULTS  FILE 

4)  EXIT 


At  this  point  the  entire  model  is  translated  with  the  creation  of  the  51  files  shown  on 
the  screen  during  execution.  Each  file  corresponds  to  a  new  time  step  and  a  new  geometric 
location  of  the  rocket.  It  is  apparent  that  running  each  file  with  PATRAN  to  obtain 
specific  data  for  a  given  node  and  component  can  become  tedious.  Therefore,  a  routine 
has  been  incorporated  into  CPI  to  allow  the  user  to  create  a  file  that  contains  only  user- 
specified  data  available  for  plotting  with  P/PLOT. 

When  the  second  option  is  chosen,  CPI  reviews  data  block  names  existing  in  the 
current  results  file.  The  user  is  asked  to  select  a  data  block  name.  CPI  classifies  data 
blocks  into  two  groups:  (1)  nodal  data,  and  (2)  element  data.  If  the  selected  data  block 
contains  nodal  data,  CPI  prompts  for  a  node  number  and  a  component  number  (1-6).  In 
general,  the  six  components  of  nodal  data  are: 

Component  1:  X-direction  translation  vector 
Component  2:  Y-direction  translation  vector 
Component  3:  Z-direction  translation  vector 
Component  4:  X-direction  rotational  vector 
Component  5:  Y-direction  rotational  vector 
Component  6:  Z-direction  rotational  vector 

If  the  selected  data  block  is  an  element  data  block,  the  user  is  prompted  for  an  element 
number  and  a  column  number. 

Next,  CPI  searches  for  the  requested  data  and,  upon  completion  of  gathering  the 
data,  requests  a  data  header  extension  of  the  current  XY-plot  data.  Up  to  38  characters  is 
acceptable.  Execution  proceeds  as  follows: 


25 


eg  A 


>»»»  ENTER  YOUR  SELECTION  <««« 


1)  WHOLE  MODEL  TRANSLATOR 

2)  XY-PLOT  DATA 

3)  ENTER  NEW  RESULTS  FILE 

4)  EXIT 


DATA  BLOCK  REVIEW 


1)  OUPV1 


ENTER  DESIRED  BLOCKNAME  NUMBER 
(OR  "0"  TO  EXIT) 

ENTER  NODE  NUMBER 

WHICH  COMPONENT  (1-6)? 

ENTER  PLOT  TITLE  EXTENSION  FOR  DATA  SET:  1 

DEFAULT  TITLE:  YDATA,  DATA  SET:  1;  NODE:  2?  COMPONENT: 


Y-TRAJECTORY  OF  ROCKET 

%%%%%  NUMBER  OF  DATA  WRITTEN  =  51  % 


0,0^,  0,0, 
'o'o'o'o 


DATA  BLOCK  REVIEW 


1)  0UPV1 


ENTER  DESIRED  BLOCKNAME  NUMBER 
(OR  "O"  TO  EXIT) 


»»»>  ENTER  YOUR  SELECTION  «««< 


1)  WHOLE  MODEL  TRANSLATOR 

2)  XY-PLOT  DATA 

3)  ENTER  NEW  RESULTS  FILE 

4)  EXIT 


CPI  Execution  completed 


At  this  point  a  file  called  XYPLOT.DAT  exists  and  upon  input  of  this  file  into  P/PLOT, 
data  is  readily  graphed  for  the  Y  coordinates  (component  2)  of  node  2  of  the  rocket  for 
each  of  the  51  time  steps  (Fig.  4). 


Fig.  4  Displacement  of  Node  2  -  Component  2 


OUTPUT  FILES 
Results  Translation 

CPI  creates  several  results  files  which  provide  a  direct  avenue  between  COS- 
MIC/NASTRAN's  analysis  algorithms  and  PATRAN's  post-processing  capabilities. 
COSMIC/NASTRAN's  input  file  (commonly  called  the  input  data  deck)  provides  the 
name  of  a  binary  file  which  CPI  interpolates.  For  example,  EXAMPLE.NAS  (the  input 
data  deck)  yields  EXAMPLE.PAT  (the  binary  file  to  be  translated  to  PATRAN- 
compatible  results  files).  Creation  of  the  binary  file  is  accomplished  by  inclusion  of  the 
appropriate  ALTER  statement  in  the  executive  control  deck  of  COSMIC/NASTRAN's 
input  data  deck.  Only  filename  extensions  are  changed  by  execution  of  COS- 
MIC/NASTRAN,  EXAMPLE.PAT  is  then  translated  by  CPI  into  appropriate  PATRAN- 
compatible  results  files.  A  description  of  each  file  follows. 

All  CPI  output  files  are  given  the  name  'Xiilnnn.EXT'  where  X  is  either  an  'S'  or  an 
'M'  depending  on  whether  the  following  two  digit  number  is  either  a  subcase  or  a  mode 
number,  respectively,  and  EXT  is  the  file  extension  named  with  respect  to  data  type. 
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XY-PIot 


Selecting  "XY-PLOT  DATA"  at  CPI's  main  menu  initiates  a  prompt  screen  which 
lists  the  various  data  blocks  found  in  COSMIC/NASTRAN's  binary  output  file.  CPI  then 
invites  the  user  to  specify  desired  data  blocks  containing  either  nodal  or  element  data  that 
are  to  be  written  to  XYPLOT.DAT,  which  is  the  input  file  for  P/PLOT.  This  file  can 
contain  one  or  more  data  sets  depending  on  how  many  times  the  user  requests  that  CPI 
write  to  this  file.  CPI  writes  only  y-data  to  XYPLOT.DAT  since  the  user  defines  an  initial 
x  and  delta-x  when  executing  P/PLOT.  A  description  of  XYPLOT.DAT's  format  is  given 
below. 

(1)  Nodal  Data  Blocks 

The  first  line  written  for  each  data  block  contains  the  default  title: 

"YDATA,  DATA  SET:  iii;  NODE  nnnn;  COMPONENT  j  -  » 

and  may  be  appended  to  allow  for  more  descriptive  titles.  This  appendage  plus  the  default 
title  must  not  exceed  80  characters.  If  more  than  80  characters  are  specified,  those  beyond 
column  80  are  truncated.  This  allows  the  user  a  suffix  of  anywhere  from  33  to  38  characters 
depending  on  the  number  of  digits  contained  in  the  default  title. 

Each  line  thereafter,  until  either  a  *Z  (end  of  file)  or  another  title  is  encountered, 
contains  nodal  data  for  components  1,  2,  3,  4,  5,  or  6,  i.e.,  X,  Y,  Z,  ex,  ©y,  or  ez, 
respectively,  as  requested  from  CPI's  prompting. 

(2)  Element  Data  Blocks 

As  with  nodal  data  blocks,  the  first  line  written  for  an  element  data  block  contains  a 
default  title  and  takes  the  form: 

"YDATA,  DATA  SET:  iii;  ELEMENT  nnnn;  COLUMN  j  -  " 

It  may  also  be  appended  as  described  above,  but  again,  the  default  title  plus  the  appendage 
must  not  exceed  80  characters  in  length  or  truncation  of  characters  beyond  column  80 
occurs.  As  with  the  nodal  data  case,  a  suffix  of  33  to  38  characters  is  allowed. 

The  remaining  lines  associated  with  the  current  data  set  contain  element  data 
available  for  plotting.  The  data  type  is  chosen  by  entry  of  a  column  number  at  CPI's 
prompt. 


SOFTWARE  AND  HARDWARE  REQUIREMENTS 

CPI  source  code  was  written  in  standard  FORTRAN  77;  however,  some  special 
features  of  VAX/VMS  FORTRAN  were  implemented.  CPI  will  run  on  any  VAX  machine 
with  a  VMS  operating  system  and  a  FORTRAN  compiler.  The  total  length  of  the  source 
code  is  approximately  1200  lines,  which  requires  104  blocks  of  disk  space. 
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ACCURACY  OF  THE  QUAD4  THICK  SHELL  ELEMENT 


William.  R.  Case,  Tiffany.  D.  Bowles,  Alicia  K.  Croft* 
and  Mark.  A.  McGinnis 
NASA/Goddard  Space  Flight  Center 


SUMMARY 


The  accuracy  of  the  relatively  new  QUAD4  thick  shell  element  is  assessed  via 
comparison  with  a  theoretical  solution  for  thick  homogeneous  and  honeycomb 
flat  simply  supported  plates  under  the  action  of  a  uniform  pressure  load.  The 
theoretical  thick  plate  solution  is  based  on  the  theory  developed  by  Reissner 
and  includes  the  effects  of  transverse  shear  flexibility  which  are  not 
included  in  the  thin  plate  solutions  based  on  Kirchoff  plate  theory.  In 
addition,  the  QUAD4  is  assessed  using  a  set  of  finite  element  test  problems 
developed  by  the  MacNeal-Schwendler  Corp.  (MSC) .  Comparison  of  the  COSMIC 
QUAD4  element  as  well  as  those  from  MSC  and  Universal  Analytics,  Inc.  (UAI) 
for  these  test  problems  is  presented.  The  current  COSMIC  QUAD4  element  is 
shown  to  have  excellent  comparison  with  both  the  theoretical  solutions  and 
also  those  from  the  two  commercial  versions  of  NASTRAN  that  it  was  compared 
to . 


INTRODUCTION 


The  QUAD4  thick  shell  element,  added  to  NASTRAN  in  1987,  is  one  of  the  most 
important  additions  to  the  program  since  the  original  writing  of  the  code. 

The  deficiencies  of  the  original  QUAD1  and  QUAD2  quadrilateral  shell  elements 
have  been  recognized  for  years  and  have  been  reported  in  the  literature.  At 
the  Goddard  Space  Flight  Center  (GSFC) ,  the  quadrilateral  shell  element  is  in 
use  in  virtaully  all  structural  analyses  of  our  spaceraft  and  related 
hardware.  Typical  applications  are  for  the  modelling  of  cylindrical  shells 
and  flat  plates  made  of  honeycomb  or  machined,  lightweighted,  metal  that  make 
up  the  structure  of  spacecraft  and  scientific  instruments.  In  some  cases 
these  models  require  that  the  effects  of  transverse  shear  flexibility  be 
included  due  to  their  thickness.  The  QUAD4  element  includes  these  effects 
and,  in  addition,  has  an  improved  isoparametric  membrane  capability  for 
in-plane  loading. 

The  purpose  of  the  study  reported  herein  is  to  assess  the  accuracy  of  the 
QUAD4  element  in  modelling  a  variety  of  situations  involving  both  solid 
cross-section  plates  as  well  as  those  constructed  of  honeycomb.  Three  goals 
of  the  study  were  to  determine: 

a)  what  is  the  rate  of  convergence  to  the  theoretical  solution  as  the 
mesh  is  refined; 
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b)  whether  the  element  exhibits  sensitivity  to  aspect  ratios 
significantly  different  than  1.0; 

c)  how  the  element  behaves  in  a  wide  variety  of  modelling 
situations,  such  as  those  included  in  the  MSC  element  test  library 
(discussed  below) . 

The  first  two  questions  were  addressed  in  the  same  manner  as  several  other 
studies  reported  by  one  of  the  authors  in  prior  NASTRAN  colloquia  (references 
1  and  2).  The  procedure  used  in  those  studies,  and  followed  here  also,  is  to 
isolate  the  effects  of  mesh  refinement  and  aspect  ratio.  That  is,  the  mesh 
refinement  study  is  done  using  elements  with  an  aspect  ratio  of  1.0.  Then., 
once  a  fine  enough  mesh  has  been  reached  such  that  the  errors  are  small,  the 
effects  of  aspect  ratio  can  be  investigated  by  keeping  the  mesh  the  same  (i.e. 
same  number  of  elements)  and  varying  the  overall  dimensions  of  the  problem, 
thus  resulting  in  each  element  aspect  ratio  changing.  Obviously,  in  order  to 
accomplish  this  latter  step  there  must  be  a  theoretical  solution  (or  some 
other  equally  acceptable  comparison  solution)  to  the  problem  with  which  to 
compare  the  finite  element  model  results.  This  is  needed  since,  at  each  step, 
a  problem  of  different  dimensions  (and  therefore  different  theoretical 
solution)  is  being  modelled. 

The  above  tests  are  important  in  that  they  show  the  rate  of  convergence  toward 
the  theoretical  solution  as  the  mesh  is  refined.  Those  tests,  however,  are 
not  sufficient  to  completely  test  the  accuracy  of  a  finite  element  since  they 
do  not  test  irregular  geometries,  or  a  variety  of  loadings  or  material 
properties.  The  MSC  has  developed  a  comprehensive  set  of  problems  for  testing 
finite  elements  in  a  variety  of  situations  (reference  3) .  The  library  of 
problems  consists  of  15  test  problems  for  the  QUAD4  element  that  cover  all  of 
the  parameters  mentioned  above.  A  test  of  the  COSMIC  QUAD4  using  these 
elements  was  reported  at  the  17th  NASTRAN  Users  Colloquium  in  1987  by  Victoria 
Tischler  of  the  Air  Force  Wright  Aeronautical  Laboratories  (AFSC)  at  Wright 
Patterson  Air  Force  Base,  Ohio,  but  was  not  included  in  the  formal 
proceedings.  Due  to  the  fact  that  it  was  not  included  in  the  formal 
proceedings,  and  also  due  to  the  fact  that  errors  in  the  QUAD4  code  for 
nonhomogeneous  plates  (to  be  discussed  later)  have  been  corrected,  the  results 
of  our  testing  of  the  latest  version  of  the  element  with  the  MSC  library  are 
include  herein. 


RESULTS  OF  MESH  AND  ASPECT  RATIO  STUDY 


For  the  mesh  and  aspect  ratio  study  a  theoretical  comparison  solution  is 
highly  desirable.  Since  the  effects  of  transverse  shear  flexibility  are 
included  in  the  QUAD4  element  formulation,  a  theoretical  solution  for 
moderately  thick  plates,  based  on  Reissner  (or  Mindlin)  thick  plate  theory  is 
also  desirable.  Such  a  solution  is  given  in  references  4  and  5  for 
rectangular  simply  supported  thick  plates  under  the  action  of  a  pressure  load. 
Thus,  this  problem  was  used  for  the  mesh  and  aspect  ratio  portions  of  the 
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Study,  Figure  defines  -tue  geometry,  coordinate  system,  boundary  conditions 
and.  loading-  fpr  the  rectangular  plate.  The  thickness  indicates  a  inoderoteiy 
thick  plate  cf  length  to  thickness  ratio  of  .29.  The  effect  of  transverse 
shear  'flexibility  is  srily  approximately  1%  on  the  maximum  displaceffiGnt  but  is 
important  in  discerning  the  qual^y  °f  the  convergence  of  the  finite  element 
results  to  the  exact  theoretical  solution.  By  exact  is  meant  the 
theoretical  basis  for  the  0UAD4  element,  'which  is  expressed  in.  the  Reissner 
thick  plate  theory .  Figure  2  shears  the  finite  element  mesh  geometry  used  in 
the  ansh  and  aspect  ratio  studies.  Due  ho  symmetry  only  one  quarter  of  the 
plate  war  modelled  The  4  x:  '4  mesh  shown  on  figure  2  is  an  example  only.;  the 
mc-sh  was  varied  during  the  .study. 

Figures  3a  -  3c  show  characteristics  of  the  theoretical  solution.  As  indicated 
in  figure  3a  Che  central  displacement  soltztion  is  represented  as  an  infinite 
series  of  hyp.erbolic  functions.  A  FORTRAN  computer  program  was  written  to 
compute  the  theoretical  solutions  for  displacements  (using  the  series  shown) 
as  well  as  stresses  (solution  not  shown).  As  m  gets  large,  where  m  is  the 
number  ef  terms  included,  in  the  series,  the  hyperbolic  functions  tend  to 
overflow  the  exponent  range  of  the  computer.  This  does  not  indicate  a  problem 
with  the  series  shewn,  as  the  hyperbolic  functions  appear  in  both  the 
numerator  and  denominator  and  their  ratio  is  numerically  stable.  However,  in 
separately  evaluating  the  numerator  and  denominator  the  overflow  problem  was 
encountered.  In  order  to  circumvent  this  problem,  the  hyperbolic  functions 
were  rewritten  in  terms  of  exponentials  allowing  the  programmed  equations,  in 
terms  o  rios  o-“  numerator  and  denominator  terms,  to  be  evaluated  without 
overflow  j/ioblems. 

Figures  3b  and  3c  show  Che  stiffness  parameters  needed  in  the  theoretical 
solution  J:oc  the  homogeneous  (i.e.  solid)  -plate  and  the  honeycomb  plate.  For 
the  honeycomb  pi Ate,  two  different  core  stiffnesses  were  investigated .  The 
scift-er  one  is  representative  of  aluminum  honeycomb  construction  that  has  been 
used  at  the  GSF-3.  The  more  flexible  one  was  chosen  because  it  represents  a 
core  flexibility  that  is  quite  low  and  was  expected  to  be  a  more  critical  test 
of  the  Qi)AD4’s  sheer  flexibility  formulation. 

The  results  of  the  mash  study,  showing  the  convergence  of  the  QUAC-4  solutions 
to  the  theoretical,  are  presented  in  tabular  form  in  tables  1-3  and  in 
graphical  form  in  figures  4-7.  Both  formats  show  %  error  in  displacement  at 
the  center  of  the  plate  as  a  function  of  mesh  refinement.  Results  are 
included  for  COSMIC  S8,  MS C  65C  and  UAI  10,0  NASTRAN .  In  the  tables  results 
for  COSMIC  version  87  i*  also  indicated  as  will  be  discussed  below.  The 
tables  merely  give  exact  numbers  (along  with  the  theoretical  displacements) 
and  the  figures  contain  the  sawe  error  information ,  but  in  graphic  form. 

Figures  4  and  5  and  table  1  are  the  results  for  the  homogeneous  plate.  The 
difference  between  the  results  in  figures  4  and  5  (and  that  in  the  two  parts 
of  table  1)  is  that  figure  4  (and  the  top  half  of  table  1)  is  lor  a  solution 
in  which  shear  flexibility  is  included  end  figure  5  (and  the  bottom  half  of 
table  1)  is  without  shear  flexibility.  These  two  situations  were  investigated 
to  rest  the  MID3  option  or,  the  PSHELL  NASTRAN  bulk  data  deck  card  which  allows 
the  effects  of  shear  flexibility  to  be  ignored  if  MID3  is  left  blank.  As  seen 
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in  figures  4  and  5  the  NASTRAN  results  converge  very  rapidly  with  mesh 
■fefineinei'X  for  COSMIC  88,  MSC  65C  and  UAI  10.0.  Table  1  contains  the  same 
information  along  with  results  for  COSMIC  87,  the  first  COSMIC  version  to 
contain  the  0UAD4  element.  As  seen,  all  versions  converge  to  less  than  0.5% 
error  for  a  mesh  size  of  8  x  8  with  the  results  without  shear  flexibility 
converging,  a  little  more  rapidly. 

Figures  6  and  7-  and  table  2  are  the  results  for  the  honeycomb  plate.  Figure  6 
(and  the  top  half  of  table  2)  are  for  the  honeycomb  plate  with  the  stiffer 
core  and  figure  7  (and  the  bottom  half  of  table  2)  are  for  the  more  flexible 
core.  As  seen  in  figures  6  and  7  the  NASTRAN  results  for  COSMIC  88  and  the 
two  commercial  NASTRAN  versions  converge  very  rapidly  for  the  two  honeycomb 
plates  as  they  did  for  the  homogenous  plate.  Table  2  contains  the  same 
information  along  with  the  results  for  COSMIC  87.  As  indicated,  the  errors  in 
the  first  version  containing  the  QUAD4  were  extremely  large  for  the  honeycomb 
plate  but,  as  reported  above,  were  quite  good  for  the  homogenous  plate.  When 
this  was  discovered  it  was  immediately  reported  to  COSMIC.  They  found  the 
problem  in  a  program  controlled  adjustable  parameter  (which,  is  used  to  avoid 
the  infamous  shear  locking  phenomena  in  earlier  thick  shell  finite  elements 
based  on  Reissner  plate  theory)  and  sent  us  a  fix  within  two  days .  After 
modifying  the  subroutine  containing  the  error,  the  results  became  that  which 
is  reported  under  the  COSMIC  88  heading  (the  same  fix  was  included  by  COSMIC 
in  the  88  release) . 

In  order  to  test  the  QUAD4's  sensitivity  to  aspect  ratio,  the  model  with  a 
12  x  12  mesh  was  run  in  which  the  plate  side  dimension  in  the  x  direction  was 
varied.  This  causes  the  element  aspect  ratio  to  vary  while  maintaining  a 
constant  mesh  in  an  attempt  to  remove  mesh  refinement  errors  from 
significantly  affecting  the  results.  As  seen  in  tables  1  and  2,  the  QUAD4 
results  with  a  12  x  12  mesh  (and  aspect  ratio  of  1.0)  have  very  little  error. 
The  results  of  the  aspect  ratio  study  are  presented  in  figures  8-10  and 
tables  3-5.  Tables  3-5  give  %  error  in  the  displacement  at  the  center  of 
the  plate  versus  aspect  ratio  for  a  model  with  a  mesh  of  12  x  12  QUAD4 
elements  (over  one  quarter  of  the  plate).  As  mentioned  above,  the  aspect 
ratio  was  varied  by  changing  the  dimension  of  the  plate  along  the  x  axis. 

Thus,  the  results  for  the  aspect  ratio  of  10  are  for  a  plate  (and  all  QUAD4 
elements)  that  is  10  times  as  long  in  the  x  direction  as  in  the  y  direction. 
Due  to  this  the  theoretical  solution  changes  with  aspect  ratio.  Figure  8  and 
table  3  are  for  the  homogenous  plate  (with  transverse  shear  flexibility)  while 
figure  9  and  table  4  are  for  the  stiff  core  honeycomb  plate  and  figure  10  and 
table  5  are  for  the  more  flexible  core  honeycomb  plate.  Investigation  of  the 
%  error  in  the  tables,  as  well  as  in  figures  8-10  show  that  the  QUAD4  has 
essentially  no  aspect  ratio  sensitivity  over  the  range  investigated. 

Based  on  the  above  results,  the  COSMIC  QUAD4  element  is  seen  to  give  very 
accurate  results  for  the  displacements  in  tl  a  problem  investigated,  both  in 
comparison  to  the  exact  theory  and  in  comparison  to  the  two  commercial 
versions  of  NASTRAN  that  we  have  at  the  GSFC.  Although  the  results  are  not 
presented  herein,  similarily  accurate  results  were  obtained  for  the  shear  and 
moment  stress  resultants  as  well. 
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RESULTS  OF  TESTING  USING  THE  MSC  ELEMENT  TEST  LIBRARY 


As  mentioned  earlier,  the  mesh  and  aspect  ratio  studies,  while  a  very  useful 
tool  in  the  evaluation  of  an  element,  do  not  test  all  of  the  important 
variables  that  affect  accuracy  in  a  finite  element  solution.  The  MSC  element 
test  library  mentioned  above  represents  a  rather  exhaustive  series  of  tests 
that  include  many  of  the  element  related  parameters  which  affect  the  accuracy 
of  a  finite  element  solution.  Reference  3  gives  a  detailed  description  of  the 
test  problems  along  with  theoretical  answers  and  the  results  of  the  testing  on 
several  MSC  elements.  The  reader  should  consult  reference  3  for  a  complete 
description  of  the  various  problems  in  the  test  series.  The  portion  of  this 
series  of  element  tests  that  relate  to  the  QUAD4  element  was  run  by  the 
authors  on  the  QUAD4  elements  contained  in  COSMIC  88,  UAI  9.8+  (not  version 
10.0  as  for  the  mesh  and  aspect  ratio  study)  and  MSC  65C.  As  the  MSC  does  in 
their  report,  the  results  are  presented  in  detail  and  also  in  a  summary  form 
in  which  the  element  is  given  a  letter  grade  of  A  through  F  based  on  the 
magnitude  of  the  error.  Table  6  shows  the  summary  results  for  the  15  tests  in 
the  series  ranging  from  a  simple  patch  test  to  modelling  of  beams  (using  the 
QUAD4  element  through  the  depth)  and  various  plates  and  shells.  The  meaning 
of  the  letter  grades  is  given  at  the  bottom  of  the  table.  As  pointed  out  in 
reference  3,  a  failing  grade  for  an  element  in  one  test  is  not  a  reason  to 
dismiss  the  element.  For  one  thing,  the  test  scores  would  improve  with  mesh 
refinement;  the  mesh  used  in  most  of  the  problems  was  quite  coarse.  Of 
importance  in  this  discussion  is  not  the  actual  grades  listed  in  table  6  but 
the  comparison  of  the  COSMIC  grades  with  those  from  the  other  two  programs. 

As  seen  in  table  6,  the  COSMIC  QUAD4  element  is  as  good  as,  or  better  than, 
those  of  the  commercial  programs.  Although  not  shown  in  table  6,  the  old 
QUAD 2  element  (included  in  reference  3)  has  a  D  or  F  grade  in  9  of  the  15 
problems.  This  is  the  reason  for  the  longstanding  need  for  an  improved  shell 
element  and  the  QUAD4  element  added  to  COSMIC  NASTRAN  clearly  fills  that  need. 
Detailed  results  for  each  of  the  problems  in  the  test  series  are  contained  in 
tables  7-12  and  are  included  for  completeness. 


CONCLUSIONS 


The  COSMIC  QUAD4  general  purpose  flat  shell  element  has  been  shown  to  be  an 
excellent  element  and  significantly  enhances  the  usefulness  of  COSMIC  NASTRAN. 
The  element  has  been  shown  to  compare  excellently  with  those  available  in  two 
commercial  versions  of  NASTRAN  that  are  currently  being  used  at  the  GSFC.  The 
addition  of  an  improved  triangular  shell  element,  anticipated  in  the  near 
future,  is  highly  desireable  as  a  companion  element  to  the  QUAD4  in  general 
analyses  of  complicated  shell  like  structures. 
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List  of  Symbols 


a,b 

t 

P 

Dj  Qi,Cs 
E 

v 

Gc 

AR,  ARe 
w 

Nx,  Ny 


=  plate  dimensions 
=  plate  thickness 
=  pressure  load 

=  plate  rigidities  (see  Figures  3b, 3c) 

=  Young’s  Modulus 
=  Poisson's  Ratio 

=  honeycomb  core  shear  modulus 
=  aspect  ratio  (ratio  of  planar  dimensions  of  plate  or  element) 

=  plate  displacement 

=  number  of  elements  in  model  of  plate  in  x,  y  directions  respectively 
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TABLE  1 
MESH STUDY 

THICK  HOMOGENEOUS  PLATE 
ELEMENT  ASPECT  RATIO  1.0 

Theoretical  Displacements 
With  Transverse  Shear  Flexibility1  3.571xl0‘5  m 

(1.406x1 O'3  in.) 

Without  Transverse  Shear  Flexibility1  3.529xl0'5  m 

(1. 390xl0-3  in.) 


%  Error 

Cosmic  Cosmic  UAI  MSC 

Mesh _ 87 _ 88  Ver.  10.0  Ver,  65C 

With  Transverse  Shear  Flexibility 


lxl 

12.03 

12.03 

12.03 

21.76 

2x2 

4.35 

4.34 

4.35 

2.54 

4x4 

1.67 

1.67 

1.67 

1.39 

8x8 

0.59 

0.60 

0.59 

0.53 

12x12 

0.39 

0.41 

0.39 

0.36 

Without  Transverse  Shear  Flexibility 

lxl 

16.90 

16.83 

16.90 

26.31 

2x2 

1.12 

1.10 

1.12 

1.67 

4x4 

0.19 

0.18 

0.19 

0.50 

8x8 

0.03 

0.00 

0.03 

0.30 

12x12 

0.00 

0.03 

0.00 

0.18 
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TABLE  2 
MESH  STUDY 

THICK  HONEYCOMB  PLATE 
ELEMENT  ASPECT  RATIO  1.0 

Theoretical  Displacements 
Gz  =  1.517xl08  N/m2 :  2.422xl0‘3  m 

(9.535xl0-2  in.) 
Gz  =  1.379xl07  N/m2 :  3.102xl0-3  m 

(1.221X10'1  in.) 


%  Error 

Cosmic  Cosmic  UAI  MSC 

Mesh  87  88  Ver.  10.0  Ver.  65C 


Gz  =  1.517xl08  N/m2  (22000.  psi) 


lxl 

747.3 

-16.31 

-7.21 

-17.98 

2x2 

589.9 

-1.17 

4.87 

3.26 

4x4 

311.4 

-0.25 

1.46 

1.19 

8x8 

103.3 

-0.06 

0.37 

0.31 

12x12 

47.9 

-0.03 

0.16 

0.14 

iz  =  1.379xl07  N/m2  (2000.  psi) 

lxl 

-6550.4 

-6.71 

10.31 

4.92 

2x2 

-5127.3 

0.26 

5.51 

4.57 

4x4 

-2689.0 

0.09 

1.42 

1.22 

8x8 

-888.5 

0.02 

0.36 

0.31 

12x12 

-412.2 

0.01 

0.16 

0.14 
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TABLE  3 

ASPECT  RATIO  STUDY 
THICK  HOMOGENEOUS  PLATE 
WITH  TRANSVERSE  SHEAR  ELEXIBILTY 
12X12  MESH 


AR 

theoretical  w, 
m  (in.) 

Cosmic 

88 

%  Error 
UAI 

Ver.  10.0 

MSC 

Ver.  65  C 

1 

3.571xl0-5 

(1.406xl0-3) 

-0.38 

0.39 

0.36 

2 

8.865xl0-5 

(3.490xl0-3) 

0.28 

0.26 

0.27 

5 

11.34x10-5 

(4.465xl0-3) 

-0.83 

-0.01 

0.05 

10 

11.38x10-5 

(4.482x10-3) 

-0.04 

-0.06 

-0.02 
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TABLE  4 

ASPECT  RATIO  STUDY 

THICK  HONEYCOMB  PLATE,  Gz=1.517  x  108  N/m2  (22000.  psi) 

12X12  MESH 


AR 

theoretical  w, 
m  (in.) 

Cosmic 

88 

%  Error 
UAI 

Ver.  10.0 

MSC 

Ver.  65C 

1 

2.422xl0'3 

(9.535xl0-2) 

0.02 

-0.16 

-0.14 

2 

5.974xl0-3 

(2.352X10'1) 

0.05 

-0.12 

-0.13 

5 

7.631x10-3 

(3.004X10'1) 

0.24 

0.13 

0.07 

10 

7.660x10-3 

(3.016X10-1) 

0.27 

0.17 

0.14 
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TABLE  5 

ASPECT  RATIO  STUDY 

THICK  HONEYCOMB  PLATE,  Gz=1.379  x  107  N/m2  (2000.  psi) 

12X12  MESH 


AR 

theoretical  w, 
m  (in.) 

Cosmic 

88 

%  Error 
UAI 

Ver.  10.0 

MSC 

Ver.  65C 

1 

3.102xl0‘3 

(1.221X10'1) 

-0.01 

-0.16 

-0.49 

2 

7.026x10-3 

(2.766X10-1) 

0.03 

-0.12 

0.23 

5 

8.785x10-3 

(3.459X10'1) 

0.20 

0.01 

0.06 

10 

8.815x10-3 
(3.470x10- !) 

0.24 

0.41 

0.14 
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TABLE  6 

SUMMARY  OF  TEST  RESULTS  FOR  QUAD4  SHELL  ELEMENTS 


Elem. 

Loading 

Test 

In 

Plane 

Out  of 
Plane 

COSMIC 

88 

UAI 

9.8+ 

MSC 

65C 

1.  Patch  Test 

X 

Irregular 

A 

A 

A 

2.  Patch  Test 

X 

Irregular 

A 

A 

A 

3.  Straight  Beam,  Extension 

X 

All 

A 

A 

A 

4.  Straight  Beam,  Bending 

X 

Regular 

B 

B 

B 

5.  Straight  Beam,  Bending 

X 

Irregular 

F 

F 

F 

6.  Straight  Beam,  Bending 

X 

Regular 

A 

A 

A 

7.  Straight  Beam,  Bending 

8.  Straight  Beam,  Twist 

X 

Irregular 

All 

A 

B 

A 

B 

B 

B 

9.  Curved  Beam 

X 

Regular 

C 

C 

C 

10.  Curved  Beam 

X 

Regular 

B 

B 

B 

11.  Twisted  Beam 

X 

X 

Regular 

A 

A 

A 

12.  Rectangular  Plate  (N=4) 

X 

Regular 

A 

A 

B 

13.  Scordelis-Lo  Roof  (N=4) 

X 

X 

Regular 

B 

B 

B 

14.  Spherical  Shell  (N=8) 

X 

X 

Regular 

A 

A 

A 

15.  Thick-Walled  Cylinder 
(nu=.4999) 

X 

Regular 

B 

B 

F 

Number  of  Failed  Tests  (D's  and  F's) 

1 

1 

2 

Grading  for  Shell  Element  Test  Results 

Grade  Requirement 

A  2%  >  Error 

B  10%  >  Error  >  2% 

C  20%  >  Error  >  10% 

D  50%  >  Error  >  20% 

F  Error  >  50% 
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TABLE  7 

PATCH  TEST  RESULTS 


Maxium  %  Error  in  Stress 


Cosmic 

UAI 

MSC 

88 

Ver  9.8+ 

Ver.  65C 

Quad4 

Quad4 

Quad  4 

Constant-Stress  Loading 

0.00 

0.00 

0.00 

Constant-Curvature  Loading 

0.00 

0.00 

0.00 
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TABLE  8 

RESULTS  FOR  STRAIGHT  CANTILEVERED  BEAM 


Normalized  Tip  Displacement* 
in  Direction  of  Loading 


Tip  Loading 

Cosmic 

UAI 

MSC 

88 

Ver.  9.8+ 

Ver.  65C 

Direction 

Quad4 

Quad4 

Quad  4 

(a)  Rectangular  Elements 

Extension 

0.996 

0.996 

0.995 

In-plane  Shear 

0.904 

0.904 

0.904 

Out-of-plane  Shear 

0.985 

0.985 

0.986 

Twist 

0.958 

0.957 

0.941 

(b)  Trapezoidal  Elements 

Extension 

1.00 

0.992 

0.996 

In-plane  Shear 

0.071 

0.071 

0.071 

Out-of-plane  Shear 

0.980 

0.979 

0.968 

Twist 

0.937 

0.934 

0.951 

(c)  Parallelogram  Elements 

Extension 

0.992 

0.992 

0.996 

In-plane  Shear 

0.080 

0.080 

0.080 

Out-of-plane  Shear 

0.986 

0.986 

0.977 

Twist 

0.895 

0.892 

0.945 

*:  Normalizing  displacement  values  listed  in  Ref.  3.  It  is  usually  a  theoretical  value. 
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TABLE  9 

RESULTS  FOR  CURVED  BEAM 


Normalized  Tip  Displacement* 
in  Direction  of  Loading 


Tip  Loading 

Cosmic 

UAI 

MSC 

88 

Ver.  9.8+ 

Ver.  65C 

Direction 

Quad4 

Quad4 

Quad  4 

In-plane  Shear 

0.834 

0.833 

0.833 

Out-of-plane  Shear 

0.971 

0.971 

0.951 

RESULTS  FOR  TWISTED  BEAM 


Normalized  Tip  Displacement* 
in  Direction  of  Loading 


Tip  Loading 

Cosmic 

UAI 

MSC 

88 

Ver.  9.8+ 

Ver.  65C 

Direction 

Quad4 

Quad4 

Quad  4 

In-plane  Shear 

0.995 

0.995 

0.993 

Out-of-plane  Shear 

0.984 

0.984 

0.985 

*:  Normalizing  displacement  values  listed  in  Ref.  3.  It  is  usually  a  theoretical  value. 
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TABLE  10 

RESULTS  FOR  RECTANGULAR  PLATE 


Normalized  Lateral  Deflection  at  Center* 


Uniform  Loac 

Concentrated  Load  f 

Cosmic 

88 

UAI 

V.  9.8+ 

MSC 

V.  65 C 

Cosmic 

88 

UAI 

V.  9.8+ 

MSC 

65C 

Nodes/ 

Edge 

Simple  Supports  1 

(a)  Aspect 

Ratio  =  1.0  j 

2 

1.01 

1.05 

0.981 

1.05 

1.04 

1.02 

4 

1.01 

1.02 

1.00 

1.02 

1.02 

1.02 

6 

1.01 

1.01 

1.00 

1.01 

1.01 

1.01 

8 

1.00 

1.01 

1.00 

1.01 

1.01 

1.01 

(b)  Aspect 

Ratio  =  5.0  | 

2 

0.986 

0.983 

1.05 

0.999 

0.989 

0.811 

4 

0.988 

0.984 

0.991 

1.02 

1.01 

0.932 

6 

0.995 

0.995 

0.997 

1.03 

1.02 

0.973 

8 

0.997 

0.997 

0.998 

1.03 

1.02 

0.989 

Clamped  Supports  | 

(a)  Aspect 

Ratio  =1.0  | 

2 

1.052 

1.046 

1.008 

0.971 

0.963 

0.994 

4 

1.038 

1.034 

1.032 

1.020 

1.015 

1.010 

6 

1.024 

1.022 

1.023 

1.027 

1.018 

1.012 

8 

1.017 

1.016 

1.016 

1.013 

1.012 

1.010 

(b)  Aspect  Ratio  =  5.0  j 

2 

1.121 

1.112 

1.314 

0.689 

0.663 

0.519 

4 

1.023 

1.019 

1.016 

0.987 

0.974 

0.863 

6 

1.013 

1.010 

1.017 

1.028 

1.019 

0.940 

8 

1.014 

1.013 

1.017 

1.034 

1.027 

0.972 

*:  Normalizing  displacement  values  listed  in  Ref.  3.  It  is  usually  a  theoretical  value. 
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TABLE  11 

RESULTS  FOR  THICK-WALLED  CYLINDER 


Normalized  Radial  Displacement* 
at  Inner  Boundary 


Poisson's 

Cosmic 

UAI 

MSC 

Ratio 

88 

Ver.  9.8+ 

Ver.  65C 

Quad4 

Quad4 

Quad4 

0.4900 

1.027 

1.027 

0.864 

0.4990 

1.032 

1.032 

0.359 

0.4999 

1.033 

1.033 

0.053 

*:  Normalizing  displacement  values  listed  in  Ref.  3.  It  is  usually  a 

theoretical  value. 

47 


TABLE  12 

RESULTS  FOR  SCORDELIS-LO  ROOF 


Normalized  Vertical  Deflection* 
at  Midpoint  of  Free  Edge 


No.  of  Spaces 

Cosmic 

UAI 

MSC 

per  Edge 

88 

Ver.  9.8+ 

Ver.  65C 

Quad4 

Quad4 

Quad4 

2 

1.450 

1.450 

1.376 

4 

1.070 

1.070 

1.050 

6 

1.030 

1.030 

1.018 

8 

1.019 

1.019 

1.008 

10 

1.015 

1.015 

1.004 

RESULTS  FOR  SPHERICAL  SHELL 


Normalized  Vertical  Deflection* 
at  Midpoint  of  Free  Edge 


No.  of  Spaces 

Cosmic 

UAI 

MSC 

per  Edge 

88 

Ver.  9.8+ 

Ver.  65C 

Quad4 

Quad4 

Quad4 

2 

1.020 

1.011 

0.972 

4 

1.043 

1.040 

1.024 

6 

1.023 

1.020 

1.013 

8 

1.010 

1.009 

1.005 

10 

1.004 

1.003 

1.001 

12 

1.000 

0.999 

0.998 

16 

0.998 

0.997 

0.997 

*:  Normalizing  displacement  values  listed  in  Ref.  3.  It  is  usually  a  theoretical  value. 
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Plate  Size:  a=1.016m(40.  in.)*  b=  1 .0 1 6m  (40.in.) 

Boundary  Conditions:  simply  supported  on  all  edges 
Loading:  pressure  load,  p=6895.  N/rrf2  ( 1.0  psi)  +Z  direction 
Thickness:  t=0.050S  m  (2.0  in.) 


*:  Variable  in  aspect  ratio  studies 


ARe=  ae/be  =  element  aspect  ratio 

Nx=  a/2ae  =  number  of  elements  in  X  direction  in  I  /4  of  plate 
N  =  b/2be  =  number  of  elements  in  Y  direction  in  1  /4  of  plate 


Fig.  3a  j 

* 


Theoretical  Solution  -  Central  Displacement 


Central  Displacement 


w(x=-y=0) 


4p  y 

3D  m»  1.3.5,  . 


[ 


+  C5  cosh(^y)  +  |iyC6  sinh(p.y) 


sin  nx 


where, 

4“-tanhM 

C6 _ ! _ 

2  cosh  am 

mrc  b  mrc 
«m=  2  a,  H=  a 
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o\  |cn 


Fig.  3b 


Theoretical  Solution  -  Homogeneous  Plate  Parameters 


Homogeneous  Plate 
Et3 

D=  12(l-v2) 


Gt, 


G- 


E 

2( 1 +v) 


E  =  6.89  x  IO'0  ,\|/m2  (10.0  x  10$  lb/in2) 
v  =  0.33 
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Fig.  3c 

Theoretical  Solution  -  Honeycomb  Plate  Parameters 


Honeycomb  Plate 


n  _  Eftf(tc+tf/?)2 

4(  1  -v2) 

Cn  =  00 
Cs  =  tcGc 

Ef  =  6.89  x  1 0 10  N/m2 
(10  x  106  lb/ in2) 

v  =  0.33 


Gc  =  1.379  x  1 07  N/m2  (2000.  Ib/in2) 
or 

1.517  x  iO8  N/m2  (22000.  Ib/in2) 
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Error  in  Displacement  at  Center  of  Plate,  % 


Fig.  4 

Error  in  Displacement  at  Center  of  Plate 

Mesh  Size  Study 

Homogeneous  Plate  with  Transverse  Shear  Flexibility 
Element  Aspect  Ratio  1.0 
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Plate,  % 


V 


Fig.  5 

Error  in  Displacement  at  Center  of  Plate 

Mesh  Size  Study 

Homogeneous  Plate  without  Transverse  Shear  Stiffness 
Element  Aspect  Ratio  1.0 


Error  in  Displacement  at  Center  of  Plate,  % 


Fig.  6 

Error  in  Displacement  at  Center  of  Plate 

Mesh  Size  Study 
Stiff  Honeycomb  Plate 
Element  Aspect  Ratio  1.0 
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Fig.  7 

Error  in  Displacement  at  Center  of  Plate 

Mesh  Size  Study 
Flexible  .Honeycomb  Plate 
Element  Aspect  Ratio  1.0 


-e-  COSMIC  88 
-b  UAIVER.  10.0 
-O  MSCVER.65C 


A  -  - 


Error  in  Displacement  at  Center  of  Plate 

Aspect  Ratio  Study 
Stiff  Honeycomb  Plate 
12x12  Mesh 


'enter  or  rlate,  % 


EIGENVALUE  COMPUTATIONS  WITH  THE  QUAD 4  CONSISTENT-MASS  MATRIX 


Thomas  A.  Butler 
Los  Alamos  National  Laboratory 


SUMMARY 


The  NASTRAN  user  has  the  option  of  using  either  a  lumped-mass 
matrix  or  a  consistent-  (coupled-)  mass  matrix  with  the  QUAD 4  shell 
finite  element.  At  the  Sixteenth  NASTRAN  Users'  Colloquium  (1988), 
Melvyn  Marcus  and  associates  of  the  David  Taylor  Research  Center 
summarized  a  study  comparing  the  results  of  the  QUAD 4  element  with 
results  of  other  NASTRAN  shell  elements  for  a  cylindrical-shell 
modal  analysis.  Results  of  this  study,  in  which  both  the  lumped- 
and  consistent-mass  matrix  formulations  were  used,  implied  that  the 
consistent-mass  matrix  yielded  poor  results.  In  an  effort  to 
further  evaluate  the  consistent-mass  matrix,  a  study  was  performed 
using  both  a  cylindrical-shell  geometry  and  a  flat-plate  geometry. 
Modal  parameters  were  extracted  for  several  modes  for  both 
geometries  leading  to  some  significant  conclusions.  First,  there 
do  not  appear  to  be  any  fundamental  errors  associated  with  the 
consistent-mass  matrix.  However,  its  accuracy  is  quite  different 
for  the  two  different  geometries  studied.  The  consistent-mass 
matrix  yields  better  results  for  the  flat-plate  geometry  and  the 
lumped-mass  matrix  seems  to  be  the  better  choice  for  cylindrical- 
shell  geometries. 


INTRODUCTION 


At  the  1988  NASTRAN  Users'  Colloquium,  results  of  a  study 
using  the  QUAD 4 ,  four-node,  shell  finite  element  for  shell 
vibrations  was  presented  (ref.  1).  This  study  indicated  that  using 
the  QUAD 4  element  with  a  consistent-mass  matrix  results  in  poor 
prediction  of  natural  frequencies  of  a  cylindrical  shell.  The 
errors  in  predicted  frequencies  were  small  for  lower  circum¬ 
ferential  harmonics  (n<4)  and  grew  from  approximately  10  per  cent 
for  the  fourth  circumferential  harmonic  to  over  20  per  cent  for  the 
sixth  circumferential  harmonic.  The  errors  seemed  to  be  relatively 
independent  of  the  longitudinal  harmonics.  The  authors  conclude 
that  the  poor  performance  is  probably  caused  by  either  a  bad 
formulation  of  the  consistent-mass  matrix  or,  more  likely,  a  coding 
error  in  the  program.  The  QUAD 4  element  is  described  in  reference 
2. 
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In  an  effort  to  determine  whether  either  of  the  above  reasons 
for  the  poor  results  is  correct,  a  study  was  undertaken  at  Los 
Alamos  to  gain  mote  insight  into  the  problem.  Earlier  studies  to 
evaluate  the  performance  of  the  element  for  static  problems 
indicate  that  the  stiffness  matrix  formulation  is  correct.  Also, 
results  reported  in  reference  1  for  the  QUAD4  element  with  a 
lumped-mass  matrix  indicate  that  the  problem  is  not  with  the 
stiffness  matrix  because  the  error  in  frequency  prediction  is  quite 
low  (less  than  4  percent  even  for  higher  circumferential  har¬ 
monics).  Of  course,  some  degradation  in  accuracy  is  expected  for 
higher  harmonics  because  the  mesh  density  can  become  a  limiting 
factor. 

As  a  first  step  in  our  evaluation,  an  independent  check  of  the 
formulation  and  coding  was  performed.  No  problems  were  found  with 
either  the  formulation  or  the  coding.  As  a  final  check,  the  mass 
matrix  for  a  single  element  oriented  at  a  skew  angle  to  the  global 
coordinate  system  was  calculated  by  hand,  and  then  results  of  the 
code  were  compared.  There  apparently  are  no  errors  in  either  the 
formulation  or  the  coding. 

A  brief  review  of  the  literature  on  the  subject  of  consistent- 
mass  matrices  does  lend  some  insight  into  the  problem.  Clough  and 
Wilson  (ref.  3)  state  that,  if  the  consistent-mass  formulation  is 
used  with  a  displacement  compatible  element,  resulting  frequencies 
are  always  larger  than  the  true  frequencies.  With  a  lumped-mass 
matrix,  the  frequencies  may  be  above  or  below  the  true  frequencies 
leading  to  the  possibility  that  use  of  the  lumped-mass  matrix 
formulation  may  result  in  more  accurate  frequency  predictions.  The 
NASTRAN  Theoretical  Manual  (ref.  4)  describes  errors  associated 
with  both  consistent-  and  lumped-mass  matrices  for  the  rod 
elements.  Fortunately,  in  this  case,  the  errors  turn  out  to  be 
opposite  in  sign  and  of  equal  magnitude  for  lower-order  terms. 
Thus,  an  accurate  mass  matrix  can  be  generated  simply  by  averaging 
the  lumped-  and  consistent-mass  matrices.  The  case  does  not  appear 
to  be  the  same  for  shell  elements.  Stavrinidis  et  al.  (ref.  5) 
propose  improved  mass  matrices  for  several  elements,  including  the 
one-dimensional  bar,  two-dimensional  membrane,  and  the  pure  bending 
beam  element.  Their  method  and  other  methods  using  so-called 
"higher-order"  mass  matrices  depend  on  the  predicted  frequencies 
being  consistently  high  or  low.  With  significant  effort,  similar 
methods  may  be  applicable  to  the  current  three-dimensional  shell 
problem.  However,  as  will  be  seen  later  in  this  paper,  solutions 
with  the  consistent-mass  matrix  for  the  QUAD4  element  can  be  either 
high  or  low,  depending  on  the  geometry  of  the  structure. 


TEST  PROBLEMS 


Two  test  problems  were  chosen  for  this  study.  The  first  was 
a  free-free  flat  plate  for  whose  natural  frequencies  we  have 
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closed-form,  analytical  solutions.  The  second  was  a  right, 
circular  cylinder.  Closed-form  solutions  do  not  exist  for  this 
geometry,  so  the  finite-difference  code  BOSOR  (ref.  6)  was  used 
with  a  fine  mesh  to  establish  the  reference  frequencies  and  mode 
shapes.  BOSOR  results  compare  favorably  with  approximate  solutions 
presented  by  Blevins  (ref.  7). 

The  flat  plate  was  a  10  by  10,  0.1-thick  square.  Its  elastic 
modulus  was  1.0  x  10s,  Poisson's  ratio  was  0.3,  and  the  density  was 
1.0.  Figure  1  shows  the  first  three  vibration  modes  of  the  plate 
with  the  theoretical  frequencies. 

The  cylindrical  shell  had  a  radius  of  300  and  a  length  of  600. 
The  material  thickness  was  3.0.  Its  elastic  modulus,  Poisson's 
ratio,  and  density  were  3.0  x  106,  0.3,  and  2.588  x  10"*.  The 
cylinder  ends  were  simply  supported  without  axial  constraint  (rigid 
diaphragm) .  Table  I  gives  the  reference  frequencies  calculated 
with  BOSOR  (ref.  6)  for  the  cylindrical  shell,  along  with  the 
approximate  solutions  given  by  Blevins  (ref.  7). 

FINITE  ELEMENT  MODELS 

Three  different  finite-element  codes  were  used  to  model  each 
of  the  two  test  problems.  The  finite-element  code  SPAR  (ref.  8) 
was  used  with  its  E43,  four-node  quadrilateral  element.  This 
element  is  based  on  a  mixed  formulation  first  proposed  by  Pian 
(ref.  9).  For  analyzing  these  problems  both  the  lumped-  and 
coupled-mass  matrices  in  the  SPAR  code  were  used.  Because  the  E43 
element  is  based  on  an  assumed-stress  function,  rather  than  an 
assumed-displacement  function,  its  coupled-mass  matrix  is  not 
"consistent."  That  is,  it  is  not  derived  from  the  same 
displacement  functions  used  in  deriving  the  stiffness  matrix.  Two 
types  of  elements  were  used  with  the  ABAQUS  finite-element  code 
(ref.  10).  The  S8R5  element  is  an  eight-node  element  that  has  only 
a  consistent-mass  matrix  option.  The  S4R5  element  is  a  four-node 
element  that  offers  only  a  lumped-mass  matrix.  Finally,  NASTRAN 
was  used  with  the  QUAD 4  element  with  both  the  lumped-  and  the 
consistent-mass  matrix  options.  In  addition,  the  problems  were 
analyzed  with  NASTRAN  using  a  matrix  that  is  the  average  of  the 
consistent-  and  lumped-mass  matrices. 

The  flat  plate  was  modeled  with  three  mesh  densities  having 
three,  five,  or  seven  nodes  along  each  edge  of  the  plate.  ABAQUS 
was  not  used  with  the  coarsest  mesh  because  that  would  have 
resulted  in  a  one-element  mesh  for  the  eight-node  S8R5.  The 
cylindrical  shell  was  also  modeled  with  three  different  mesh 
densities.  These  meshes  had  5,  9,  or  17  nodes  on  each  edge.  For 
the  ABAQUS  eight-node  element,  fewer  total  nodes  were  present 
because  of  the  lack  of  the  middle  node.  For  this  study,  only  one 
eighth  of  the  shell  was  modeled,  and  symmetry  conditions  were  used 
on  all  boundaries.  Thus,  only  the  even  circumferential  and  odd 
longitudinal  harmonics  were  determined. 
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All  the  NASTRAN  solutions  were  obtained  using  the  FEER 
eigenvalue  extraction  method.  The  mass-orthogonality  test 
parameter  was  0.0001  for  the  analyses. 


RESULTS 


Results  for  the  flat  plate  are  shown  in  figures  2  through  4. 
In  these  figures,  the  horizontal  axes  show  the  number  of  nodes  on 
each  side  of  the  square  mesh  and  the  vertical  axis  is  the  natural 
log  of  the  absolute  value  of  the  error  in  predicted  frequency.  The 
error  is  simply  the  ratio  of  the  calculated  frequency  to  the 
theoretical  frequency  minus  1.0.  For  reference,  a  plotted 
In  (error)  of  -4.6  is  approximately  1.0  per  cent  in  error  in 
absolute  frequency  determination.  A  plotted  value  of  -8.0  is 
roughly  equivalent  to  0.03  per  cent  error. 

The  data  points  labeled  ’'lumped,"  "consistent,"  and  "average" 
are  all  for  the  NASTRAN  QUAD4  element.  Study  of  the  results 
reveals  some  definite  patterns.  As  might  be  expected,  the 
consistent-mass  matrix  always  outperforms  the  lumped-mass  matrix. 
However,  the  rates  of  convergence  seem  to  be  approximately  the 
same.  The  SPAR  results  that  were  obtained  by  using  the  coupled- 
mass  matrix  are  consistently  better  than  the  NASTRAN  results. 
However,  the  convergence  pattern  is  not  smooth  and,  for  all  cases, 
the  SPAR  E43  element  with  its  coupled-mass  matrix  yields  better 
answers  with  the  intermediate,  rather  than  the  fine,  mesh.  This 
result  is  somewhat  disturbing,  although,  in  all  cases,  the  errors 
were  small.  The  ABAQUS  S8R5  element  also  gives  slightly  better 
results  than  does  the  NASTRAN  QUAD 4  element.  For  the  flat  plate, 
the  elements  with  the  consistent-mass  matrix  formulation  always 
overpredicted  the  frequencies  and  those  with  the  lumped-mass 
matrices  always  underpredicted  the  frequencies. 

Results  for  the  cylinder  are  not  as  clear  as  for  the  flat 
plate.  Figures  5  through  7  show  the  frequency-convergence 
characteristics  for  the  elements  that  are  being  considered  for 
three  different  modes.  These  involve  the  second,  fourth,  and  sixth 
circumferential  harmonics  (n=2,4,and  6)  and  the  first  longitudinal 
harmonic  (m=l).  The  most  striking  observation  is  that,  for  the 
QUAD 4  element,  the  lumped-mass  matrix  is  now  outperforming  the 
consistent-mass  matrix.  This  observation  seems  to  confirm  the 
result  of  Marcus  (ref.  1).  To  illustrate  the  point,  data  from 
reference  1  have  been  added  to  the  figures.  Here,  the  definition 
of  the  ordinate  axis  has  to  be  qualified.  In  reference  1,  a  13 
node  by  37-node  mesh  was  used  in  modeling  one  half  of  a  cylinder. 
This  becomes  a  7-node  by  19-node  mesh  when  an  eighth  of  the 
cylinder  is  considered,  as  is  the  case  for  this  study.  Because, 
for  the  modes  presented,  only  the  first  longitudinal  harmonic  is 
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present,  we  can  loosely  define  this  as  a  7  by  7  mesh  and  plot  it 
as  such  on  our  figures. 

Only  for  the  lower,  second  harmonic  (fig.  5)  does  the  ABAQUS 
S8R5  element  outperform  the  NASTRAN  lumped-mass  element.  For  all 
three  modes,  the  QUAD4  with  the  lumped-mass  matrix  yields  the  best 
results.  Its  deviation  from  the  QUAD4  with  the  consistent-mass 
matrix  increases  with  higher  circumferential  harmonics.  The  SPAR, 
E43  element  with  a  coupled-mass  matrix  tends  to  follow  the  QUAD4 
element  closely  for  these  modes. 

Except  for  a  few  cases,  the  frequencies  were  overpredicted 
for  consistent-,  lumped-,  and  coupled-mass  matrices  for  the 
cylindrical-shell  problem.  The  exceptions  were  the  QUAD4  lumped- 
mass  matrix  and  the  E43  coupled-mass  matrix  for  the  finest  mesh  for 
the  second  and  third  modes  considered  here. 

Frequency  is  not  the  only  parameter  that  should  be  considered 
for  modal  analyses.  The  other  is,  of  course,  the  mode  shape.  One 
method  of  comparing  mode  shapes  is  to  compare  calculated 
generalized  masses  for  the  solutions  using  the  different  elements 
being  considered.  Another  is  to  use  a  parameter  frequently 
calculated  when  comparing  calculated  mode  shapes  with 
experimentally  measured  mode  shapes.  This  parameter  is  called  the 
mode  shape  correlation  coefficient  (MSCC)  and  is  described  in 
reference  11.  It  essentially  provides  a  measure  of  the  least- 
squares  deviation  of  the  points  being  compared  from  a  straight-line 
correlation.  Both  these  measures  were  used  in  comparing  solutions 
for  the  n=8,  m=5  mode  for  the  cylinder  being  considered  here. 
Results  of  these  comparisons  are  summarized  in  table  II,  along  with 
comparisons  of  the  predicted  frequencies.  The  predicted 
frequencies  are  normalized  using  the  BOSOR  code  results  as  the 
baseline.  The  generalized  masses  were  normalized  using  the 
theoretical  value  obtained  by  direct  integration  of  the  square  of 
the  analytically  perfect  mode  shape  multiplied  by  the  material 
density.  For  the  MSCC  comparisons,  mode  shapes  predicted  by  BOSOR 
were  used  as  the  "correct"  shape. 

A  study  of  the  results  summarized  in  table  II  shows  again  that 
the  lumped-mass  matrix  provides  better  frequency  predictions  than 
does  the  consistent-mass  matrix  for  the  NASTRAN  QUAD 4  shell 
element.  Note  that  for  the  9-node  by  9-node  mesh,  the  error  for 
the  consistent-mass  matrix  is  over  30  per  cent.  A  finer  mesh  (17 
by  17)  with  the  consistent-mass  matrix  provides  better  frequency 
approximations,  but  the  prediction  is  still  not  as  good  as  for  the 
lumped-mass  matrix  with  a  coarser  mesh.  The  generalized  mass  is 
in  considerable  error  for  both  QUAD 4  cases  in  which  the  consistent- 
mass  matrix  is  used. 

The  generalized  mass  is  a  much  more  sensitive  measure  of  mode 
shape  error  than  the  MSCC,  as  evidenced  by  data  for  the  ABAQUS 
results  that  used  the  S8R5  element.  Here  the  MSCC  is  quite  close 
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to  1.0  for  the  9-node  by  9-node  mesh,  but  the  generalized  mass  is 
over  30  per  cent  in  error.  As  the  mesh  is  refined,  the  generalized 
mass  improves,  but  it  is  still  not  as  accurate  as  for  the  QUAD 4 
when  a  lumped-mass  matrix  is  used.  Note  that  the  performance  of 
the  ABAQUS  S4R5  element  compares  favorably  with  the  NASTRAN  QUAD4 
element. 

The  SPAR  E43  element,  which  performed  nearly  as  well  as  the 
QUAD4  in  predicting  frequencies  for  all  the  shell  modes  considered 
in  this  study,  apparently  predicts  both  the  frequency  and 
generalized  mass  accurately  if  the  coupled-mass  matrix  is  used. 
However,  somewhat  unexpectedly,  this  element  does  not  perform  quite 
so  well  with  a  lumped-mass  matrix.  In  this  case,  the  frequency  is 
predicted  accurately  but  the  mode  shape  has  considerable  error 
associated  with  it,  as  evidenced  by  the  underprediction  of  the 
generalized  mass. 


CONCLUSIONS 

Among  the  elements  considered  in  this  study,  the  NASTRAN  QUAD4 
element  with  a  lumped-mass  matrix  seems  to  be  the  best  choice  when 
the  geometry  is  a  cylindrical  shell.  A  general  rule  seems  to  be 
that,  for  any  element  considered  here,  consistent-mass  matrices 
should  be  avoided  for  this  particular  geometry.  On  the  other  hand, 
for  flat-plate  geometries,  the  consistent-mass  matrix  outperforms 
the  lumped-mass  matrix. 

These  conclusions  imply  that  choices  are  difficult  when 
modeling  geometries  that  deviate  from  the  simple  geometries 
considered  here.  It  is  possible  that  an  alternate  method  of 
deriving  the  mass  matrix,  such  as  the  SPAR  coupled-mass  matrix, 
would  generate  a  result  that  would  be  more  generally  applicable. 
Note  that  it  seems  to  perform  well  for  both  geometries.  However, 
for  the  present,  if  the  geometry  is  predominantly  cylindrical,  the 
lumped-mass  matrix  should  always  be  used  with  the  NASTRAN  QUAD 4 
element. 
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TABLE  I 

PREDICTED  FREQUENCIES  (Hz)  FOR  CYLINDRICAL  SHELL  USING  BOSOR  (BLEVINS) 

FCR  EVEN  CIRCUMFERENTIAL  HARMONICS  (n)  AND  ODD  LONGITUDINAL  HARMONICS  (m). 


\m 

n\ 

1 

3 

5 

2 

19.61  (21.82) 

47.97  (48.61) 

54.69  (54.83) 

4 

7.92  (8.27) 

33.27  (33.85) 

47.05  (47.30) 

6 

7.32  (7.59) 

23.64  (24.00) 

39.56  (39.83) 

8 

10.76  (11.68) 

20.63  (20.94) 

35.18  (35.47) 

TABLE  II 


COMPARISON  OF  FREQUENCY,  GENERALIZED  MASS,  AND  MODE  SHAPE  PREDICTED 
BY  VARIOUS  FINITE-ELEMENT  MODELS  FOR  CYLINDRICAL-SHELL  MODE  n=8,  m=5. 


Nodes/side 

Computer  code/ 

element 

Mass 

matrix 

Normalized 

frequency 

Mode  shape 

correlation  coef. 

Normalized 

generalized  mass 

9 

SPAR/E43 

coupled 

1.032 

0.9995 

1.011 

9 

SPAR/E43 

lumped 

1.046 

0.9853 

0.810 

9  ABAQUS/S8R5  consistent  0.982 _ 0.9998 _ 0.677 

17  ABAQUS/S8R5  consistent  1.010  0.971 


9  NASTRAN/QUAD4  lumped  1.014  0.9991  1.010 


9  NASTRAN/QUAD4  consistent  1.336 _ 0.603 

17  NASTRAN/QUAD4  consistent  1.066  0.876 


mode  1  mode  2  mode  3 

f= 0.2055  f=0.3015  f=0.3722 


Figure  1.  Mode  shapes  and  frequencies  (Hz)  for  flat  plate. 


lumped  consistent  average  SPAR  (coupled) 

— B—  —A—  •■•©■■■ 


Figure  2.  Frequency  convergence  for  mode  1  of  a  flat  plate 
as  a  function  of  mesh  density  for  different  finite  elements. 
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ABAQUS 

(consistent)s^^ 


4  6  8  10  12  14  16  18 

nodes  per  side 

lumped  consistent  average  SPAR  (coupled) 

— S—  —A-  O 

Figure  3.  Frequency  convergence  for  inode  2  of  a  flat  plate 
as  a  function  of  mesh  density  for  different  finite  elements. 


nodes  per  edge 

lumped  consistent  average  SPAR  (coupled) 
- 0 -  -••A--’  .  Q..  - - 


Figure  4.  Frequency  convergence  for  mode  3  of  a  flat  plate 
as  a  function  of  mesh  density  for  different  finite  elements. 


lumped  consistent  average  SPAR  (coupled) 

— S —  — & —  --O'  — — 

Figure  5.  Frequency  convergence  for  mode  n=2,  m=l  of  cylindrical 
shell  as  a  function  of  mesh  density  for  different  finite  elements. 


nodes  per  side 

lumped  consistent  average  SPAR  (coupled) 
— S—  —A—  O 


Figure  6.  Frequency  convergence  for  mode  n=4,  m=l  of  a  cylindrical 
shell  as  a  function  of  mesh  density  for  different  finite  elements. 
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lumped  consistent  average  SPAR  (coupled) 


Figure  7.  Frequency  convergence  for  mode  n=6,  n=l  of  a  cylindrical 
shell  as  a  function  of  mesh  density  for  different  finite  elements. 
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ABSTRACT 


OOSMIC/NASTRAN,  as  it  is  supported  and  maintained  by  COSMIC,  runs  on  four 
main-frame  computers  -  CDC,  VAX,  IEM  and  UNIVAC.  OOSMIC/NASTRAN  on  other 
computers,  such  as  CRAY,  AMDAHL,  PRIME,  CONVEX,  etc. ,  is  available  commercially 
from  a  number  of  third  party  organizations.  All  these  computers,  with  their  own 
one-of-a-kind  operating  systems,  make  NASTRAN  machine  dependent.  The  job  control 
language  (JCL) ,  the  file  management,  and  the  program  execution  procedure  of  these 
computers  are  vastly  different,  although  95  percent  of  NASTRAN  source  code  was 
written  in  standard  ANSI  FORTRAN  77. 

The  advantage  of  the  UNIX  operating  system  is  that  it  has  no  machine 
boundary.  UNIX  is  becoming  widely  used  in  many  workstations,  mini's,  super-PC's, 
and  even  some  main-frame  computers.  NASTRAN  for  the  UNIX  operating  system  is 
definitely  the  way  to  go  in  the  future,  and  makes  NASTRAN  available  to  a  host  of 
computers,  big  and  small. 

Since  1985,  many  NASTRAN  improvements  and  enhancements  were  made  to  conform 
to  the  ANSI  FORTRAN  77  standards.  A  major  UNIX  migration  effort  was  incorporated 
into  COSMIC  NASTRAN  1990  release.  As  a  pioneer  work  for  the  UNIX  environment,  a 
version  of  COSMIC  89  NASTRAN  was  officially  released  in  October  1989  for  DEC 
UITRIX  VAXstation  3100  (with  VMS  extensions) .  A  COSMIC  90  NASTRAN  version  for  DEC 
ULTRIX  DECstation  3100  (with  RISC)  is  planned  for  April  1990  release.  Both 
workstations  are  UNIX  based  computers.  The  COSMIC  90  NASTRAN  will  be  made 
available  on  a  TK50  tape  for  the  DEC  ULTRIX  workstations.  Previously  in  1988,  an 
88  NASTRAN  version  was  tested  successfully  on  a  SiliconGraphics  workstation. 


INTRODUCTION 


The  advantage  of  AT&T's  UNIX  operating  system  is  that  it  is  an  "open 
system",  hardware  independent,  single  and  multiuser  system,  powerful,  versatile, 
and  reliable.  This  "open  system",  which  may  appear  under  different  names  such  as 
ULTRIX,  XENIX,  SunOS,  AIX  etc.,  is  becoming  the  standard  software  today  for  the 
fast-growing  market  of  workstation  computers.  Even  IEM  is  going  to  adopt  UNIX  for 
its  forthcoming  workstations.  As  many  more  computers  are  designed  to  run  under  the 
UNIX  banner,  these  newcomers  are  getting  cheaper,  faster,  and  more  powerful.  The 
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result:  unprecedented  price  competition  that's  making  UNIX  another  word  for  cheap 
computing.  The  migration  of  NASTRAN  to  the  UNIX  "open  system",  is  definitely  the 
way  to  go. 


THE  EARLY  DEVEIOMENT 


NASTRAN  is  written  mainly  in  FORTRAN  language.  Only  about  five  percent  of 
the  source  codes  are  machine-dependent.  The  early  stage  of  migration,  started  in 
1984-1985,  was  a  move  towards  ANSI  FORTRAN  77,  which  is  a  standard  FORTRAN 
compiler  for  all  UNIX  based  computers.  In  this  early  stage  of  development,  the 
NASTRAN  UNIVAC  version  was  moved  from  the  'FOR'  compiler  to  the  'FIN'  compiler, 
and  the  CDC  version  from  FORTRAN  4  to  FORTRAN  5.  The  VAX  NASTRAN  had  been 
maintained  as  a  separate  version  until  the  1984  release.  This  release  shared  the 
machine  independent  source  code  with  the  other  computers  (IBM,  CDC  and  UNIVAC) . 


TEST  ON  SiliconGraphics  WORKSTATION 


A  NASTRAN  test  program,  based  on  COSMIC/NASTRAN  88  VAX  release,  was 
converted  and  ran  "successfully"  on  a  SiliconGraphics  workstation.  Only 
occasionally  this  test  program  failed  in  some  NASTRAN  dynamic  problems.  Several 
UNIX  job  control  languages  (JCXs)  were  written  to  compile,  link  edit,  and  execute 
NASTRAN  for  this  test  program.  These  JCIs  played  an  important  part  in  the  success 
of  the  SiliconGraphics  pilot  NASTRAN  test.  With  further  refinement  and  improvement 
(done  in  1989) ,  the  JOs,  applicable  to  all  UNIX  based  computers,  play  an 
important  role  in  the  NASTRAN  migration  to  UNIX.  The  JCL  to  execute  a  NASTRAN  job 
(cold  start,  restart  or  substructuring)  is  indeed  very  user  friendly. 

This  SiliconGraphics  test  program  was  also  used  to  identify  and  verify 
efficiency  improvements  of  the  NASTRAN  source  code.  The  UNIX  utility  profiler, 
prof,  was  used  for  timing  studies  of  the  codes  needing  efficiency  improvement. 
These  studies  resulted  in  over  30  percent  speed  improvement  of  the  VAX  NASTRAN 
version.  The  other  NASTRAN  versions  were  also  benefited. 

All  changes  that  were  required  to  make  this  SiliconGraphics  test  program 
successful,  were  incorporated  into  the  machine  independent  NASTRAN  source  codes. 


VAX  NASTRAN 


The  VAX  version  of  NASTRAN  is  written  entirely  in  FORTRAN  language. 
Hardware-wise,  VAX  and  many  UNIX  based  computers  are  quite  similar.  They  are 
virtual  memory  computers  with  32-bit  word  architecture.  The  file  management 
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systems  are  quite  similar.  The  VAX  FORTRAN  version  is  the  natural  choice  for  the 
migration  of  NASTRAN  to  the  UNIX  system. 

The  NASTRAN  GINO  (General  Input  aNd  Output  file  management)  package  of  the 
VAX  version  has  gone  through  extensive  revision  and  improvement  in  1987-1988.  Many 
I/O  processes  have  been  shortened  and  streamlined.  The  packing  and  unpacking  of 
matrix  data  were  improved  and  speeded  up.  The  UNIX  based  computers  have  therefore 
benefited  from  previous  VAX  improvements.  (The  improved  VAX  GINO  and  matrix 
packing/unpacking  was  also  tested  successfully  on  an  IBM  3084  machine) 

The  VAX  89  NASTRAN  release  was  compiled  and  linked  successfully  on  a  DEC 
UITRIX  VAXstation  3100  (with  VMS  extensions) ,  using  the  UNIX  JCIs  from  the 
SiliconGraphics  test  program.  Only'  one  subroutine,  CKUTM  that  obtains  the  CHJ 
time  from  the  computer  system,  needed  modification.  All  119  NASTRAN  demonstration 
problems,  plus  20  more  user  problems,  ran  successfully.  This  NASTRAN  UNIX  version 
was  officially  released  on  a  TK50  tape  in  October  1989. 

The  DEC  UITRIX  DECstation  3100  (with  RISC,  Reduced  Instruction  Set  Chip) 
required  additional  modification  of  the  NASTRAN  source  codes.  (See  next 
paragraph.)  Occasionally  this  version  failed  in  some  dynamic  problems,  exhibiting 
the  same  symptom  as  that  of  the  SiliconGraphics  test  program.  There  will  be  no 
official  89  release  of  this  UNIX  version.  Presently,  it  is  planned  to  have  a  90 
NASTRAN  release  for  UNIX  based  computers  with  RISC  processors. 


UNIX/FORTRAN  REQUIREMENTS 


The  ANSI  FORTRAN  77  is  the  standard  FORTRAN  compiler  for  all  UNIX  based 
computers  and  workstations.  However,  small  differences  may  exist  among  ANSI 
FORTRAN  77  compilers  from  different  manufacturers.  The  1990  OOSMIC/NASTRAN 
incorporates  many  known  specifications  that  are  required  by  various  ANSI  FORTRAN 
77  compilers.  The  changes  involve: 

a.  External  declaration  of  bit-shifting  functions  (LSHIFr  and  RSHIFT) , 
the  logical  functions  (ANDF  and  ORF) ,  and  complement  function 
(OOMPLF) ,  to  avoid  system  functions  of  the  same  names. 

b.  Standardization  of  OPEN,  READ  and  WRITE  commands  for  direct-access 
files. 

c.  Removal  of  octal  and  hexadecimal  constants  from  FORTRAN  executable 
source  code. 

d.  Elimination  of  jumping  into  an  inner  do-loop,  which  was  previously 
allowed  via  an  ASSIGN  statement. 

e.  Dimension  of  one  for  all  open  core  arrays. 

f.  Alignment  of  all  open  core  arrays. 
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The  last  change  (f)  above  is  the  most  tricky  process  for  UNIX  based 
computers,  particularly  those  with  RISC.  Throughout  the  NASTRAN  source  code, 
several  hundred  labelled  commons  are  used  for  the  open  core  space.  (NASTRAN  has  no 
dimension  limit.  The  open-ended  open  core  is  used  as  scratch  space  for  internal 
computations  and  storage) .  The  other  mainframe  computers,  particularly  the  non¬ 
virtual  CDC  and  UNIVAC  machines,  require  this  open  core  space  to  have  a  unique 
name  in  a  subroutine  or  group  of  subroutines,  such  that  the  open  core  space  can  be 
positioned  strategically  in  the  executable  overlay  program.  To  cortpramise  among 
the  virtual  (UNIX  based  computers,  VAX  and  I  EM)  and  the  non-virtual  memory 
computers  (that  require  program  overlays),  a  block  data  routine,  ZZCORE.f,  was 
written  to  be  used  only  for  the  UNIX  based  computers.  All  the  open  core  labelled 
commons  are  included  in  this  block  data  routine.  The  labelled  common  /XNSTRN/  must 
be  the  very  first  in  the  list,  and  /ZZZZZZ/  must  be  the  very  last.  These  first  and 
last  labelled  common  requirements  must  be  true  not  only  in  the  FORTRAN  source 
code,  but  also  in  the  compiled  relocatable  (or  object)  program.  The  user  could  use 
the  UNIX  command  'nm  -n  ZZOORE.o'  to  verify  that  /XNSTRN/  and  /ZZZZZZ/  are 
positioned  correctly.  If  they  are  not,  something  must  be  done  to  get  the  NASTRAN 
open  core  alignment  correct.  It  is  for  this  reason  (too  many  labelled  commons  in 
one  subroutine)  that  the  DEC  UITRIX  DECstation  3100  (RISC)  uses  two  block  data 
routines:  ZZCORE.f  for  NASTRAN  links  1  through  14,  and  ZZKORE.f  for  NASTRAN  link 
15.  It  is  also  for  this  reason  that  a  C-program,  SOROBJ.c,  is  included  in  the 
UNIX  NASTRAN  release  tape,  to  sort  the  open  core  labelled  commons  in  ZZCORE.o  (a 
relocatable  file) ,  only  if  all  other  efforts  fail  to  obtain  the  proper  alignment. 


CONCLUSION 


The  90  OOSMIC/NASTRAN  release  incorporates  many  changes  as  required  by  the 
UNIX  based  computers  and  workstations.  With  a  set  of  proven  user  friendly  UNIX 
JCDs,  it  should  run  successfully  on  many  UNIX  based  computers  presently  available, 
or  still  on  the  vendors'  drawing  boards.  (FORTRAN  compile  and  link  edit  are 
required.)  Of  course,  this  90  OOSMIC/NASTRAN  release  will  continue  to  operate  as 
before  on  the  I  EM,  VAX,  CDC  and  UNIVAC  mainframes.  This  is  a  version  that  bridges 
from  the  old  world  of  proprietary  and  limited  operating  systems  to  the  new  UNIX 
world  of  "open  system". 
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OBTAINING  EIGENSOLUTIONS  FOR  MULTIPLE  FREQUENCY 
RANGES  IN  A  SINGLE  NASTRAN  EXECUTION 


by 

P.  R.  Pamidi 
RPK  Corporation 
Columbia,  Maryland 

and 


W.  K.  Brown 
RPK  Corporation 
Hayes,  Virginia 


SUMMARY 

A  novel  and  general  procedure  for  obtaining  eigenvalues  and  eigenvectors  for  multiple  fre¬ 
quency  ranges  in  a  single  NASTRAN  execution  is  presented.  The  scheme  is  applicable  to  normal 
modes  analyses  employing  the  FEER  and  Inverse  Power  methods  of  eigenvalue  extraction.  The 
procedure  is  illustrated  by  examples. 


INTRODUCTION 

NASTRAN  currently  offers  four  methods  for  real  eigenvalue  extraction.  They  are  the  Tridi¬ 
agonal  or  Givens  method,  the  Tridiagonal  Reduction  or  FEER  method,  the  Inverse  Power  method 
and  the  Determinant  method  (see  Section  10  of  Reference  1  for  details). 

The  Givens  method  is  a  transformation  method  that  computes  all  of  the  eigenvalues  in  a  prob¬ 
lem,  in  addition,  eigenvectors  corresponding  to  a  specified  range  of  frequencies  or  to  a  specified 
number  of  lowest  eigenvalues  can  also  be  computed.  The  FEER  method  is  also  a  transformation 
method  that  allows  for  the  extraction  of  a  specified  number  of  eigensolutions.  It  requires  the 
specification  of  a  “shift  point”  or  frequency  around  which  the  eigensolutions  are  desired.  The 
Inverse  Power  and  Determinant  methods  are  both  tracking  methods  that  allow  for  the  extraction  of 
a  specified  number  of  eigensolutions.  They  both  require  the  specification  of  a  frequency  range  for 
which  the  eigensolutions  are  desired. 

When  all,  or  almost  all,  of  the  eigenvalues  in  a  problem  are  required,  the  Givens  method  is  the 
generally  the  most  efficient  method  to  use  because  the  total  effort  is  not  highly  dependent  upon  the 
number  of  eigenvalues  that  are  extracted.  However,  when  the  order  of  the  problem  exceeds  a  few 
hundred,  this  method  may  require  prohibitively  time-consuming  out-of-core  operations,  thereby 
losing  its  efficiency. 
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If  an  user  is  not  interested  in  obtaining  all  of  the  modes  in  a  problem,  but  only  in  a  smaller 
number  of  eigensolutions  around  a  certain  frequency  or  within  a  certain  frequency  range,  the  FEER, 
the  Inverse  Power  and  the  Determinant  methods  are  the  obvious  choices.  The  computations  in  all 
of  these  methods  are  proportional  to  the  number  of  eigensolutions  extracted.  The  FEER  method  is 
probably  the  most  efficient  of  the  three.  It  is  quite  effective  in  obtaining  eigensolutions  around  the 
selected  frequency.  It  is  also  very  efficient  even  when  out-of-core  operations  are  involved.  The 
results  obtained  by  the  Inverse  Power  and  Determinant  methods,  on  the  other  hand,  are  very 
susceptible  to  the  number  of  estimated  roots  specified  for  a  given  frequency  range  (field  5  on  the 
EIGR  bulk  data  card;  see  Reference  2).  When  the  specified  number  for  the  estimated  roots  is  larger 
than  the  actual  number  of  roots  in  that  range,  these  methods  are  apt  to  yield  many  lower  frequencies 
outside  the  specified  frequency  range.  The  Determinant  method  is  the  least  efficient  of  all  of  the 
methods  and  will,  therefore,  not  be  considered  any  further  for  the  purpose  of  this  paper. 


CURRENT  PROCEDURE  FOR  OBTAINING  EIGENSOLUTIONS 
FOR  MULTIPLE  FREQUENCY  RANGES 

There  are  practical  situations  in  which  an  user  may  be  interested  in  obtaining  eigensolutions  for 
multiple  frequency  ranges,  w  ith  one  frequency  range  quite  distinct  and  apart  from  another  frequency 
range.  Complex  configurations  involving  control  systems,  experimental  setups  and  structural 
subsystems  are  examples  of  such  situations. 

None  of  the  eigenvalue  extraction  methods  discussed  earlier  can  accomplish  the  above  objec¬ 
tive  directly.  So,  if  the  user  wishes  to  obtain  eigenvalues  for  more  than  one  range  of  frequencies  in 
such  cases  as  the  above,  he  can  accomplish  it  at  present  in  one  of  two  ways.  The  first  way  is  for  the 
user  to  make  a  single  NASTRAN  execution  with  a  large  frequency  range  (or  a  shift  point  in 
conjunction  with  a  large  number  of  desired  roots)  so  as  to  encompass  all  of  the  frequencies  in  the 
ranges  of  interest.  However,  this  will  not  be  very  cost  effective  if  the  frequency  ranges  of  interest 
are  w  idely  separated.  The  alternative  way  is  for  the  user  to  perform  multiple  NASTRAN  executions, 
one  execution  for  each  range  of  frequencies,  effectively  utilizing  the  APPEND  feature  (see  Section 
9.2.2  in  Reference  1  and  Section  2.3.7  in  Volume  2,  Reference  2).  However,  this  latter  procedure 
involves  checkpoint/restart  runs  and  is  rather  cumbersome  for  the  purpose. 


PROPOSED  PROCEDURE  FOR  OBTAINING  EIGENSOLUTIONS 
FOR  MULTIPLE  FREQUENCY  RANGES 

The  above  objective  of  obtaining  eigensolutions  for  multiple  frequency  ranges  can  be  accom¬ 
plished  in  a  single  NASTRAN  execution  by  an  innovative  procedure  that  involves  the  use  of  DMAP 
ALTERS  in  conjunction  with  certain  specific  input  data  requirements.  This  procedure  involves 
performing  a  normal  modes  analysis  employing  multiple  subcases  and  using  the  FEER  method  or 
the  Inverse  Power  method  (whichever  is  preferred).  Each  subcase  is  setup  so  as  to  obtain 
eigensolutions  in  a  specified  frequency  range.  The  Final  results  of  the  analysis  will  contain  the 
eigensolutions  obtained  in  all  of  the  specified  frequency  ranges.  This  procedure  is,  in  essence,  a 
novel  application  of  the  APPEND  feature  referred  to  above.  The  important  difference  is  that,  while 
the  APPEND  feature  was  originally  conceived  to  be  employed  in  a  checkpoint/restart  env  ironment, 
the  proposed  procedure  accomplishes  this  in  a  single  NASTRAN  execution. 
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The  DMAP  ALTER  package  required  for  the  above  procedure  is  given  in  Appendix  A.  The 
details  of  the  input  data  requirements  and  the  output  obtained  from  the  analysis  are  discussed  below. 

Executive  Control  Deck 


The  user  should  employ  the  DMAP  ALTER  package  given  in  Appendix  A  either  by  explicitly 
including  it  or  by  referencing  it  via  a  READFILE  card  in  the  Executive  Control  Deck.  Note  that, 
for  every  additional  subcase  beyond  the  first  subcase,  the  ALTER  package  involves  a  pair  of  OFP 
and  READ  modules. 

Case  Control  Deck 


The  user  should  have  as  many  subcases  in  the  Case  Control  Deck  as  the  number  of  distinct 
frequency  ranges  for  which  he  wishes  to  obtain  eigensolutions.  Each  subcase  must  have  a  separate 
METHOD  request.  The  METHOD  request  in  each  subcase  references  a  distinct  EIGR  card  in  the 
Bulk  Data  Deck  that  either  implies  (in  the  case  of  the  FEER  method)  or  defines  (in  the  case  of  the 
Inverse  Power  method)  a  distinct  frequency  range. 

All  output  requests  and  constraint  specifications  must  be  above  the  subcase  level.  Thus,  the  only 
difference  between  one  subcase  and  another  must  be  the  different  METHOD  that  they  request.  Also, 
since  the  final  results  include  eigensolutions  for  all  of  the  subcases,  PLOT  requests  should  not  make 
explicit  references  to  subcase  numbers. 

Bulk  Data  Deck 


In  addition  to  the  required  modeling  (geometry,  constraints,  etc.)  data,  the  Bulk  Data  Deck 
should  have  as  many  EIGR  cards  as  the  number  of  subcases  employed  (and  the  corresponding 
METHOD  requests)  in  the  Case  Control  setup. 

When  using  the  FEER  method,  the  EIGR  bulk  data  card  for  each  METHOD  request  requires 
the  specification  of  a  shift  point  or  frequency  that  indicates  the  center  of  a  frequency  range  as  well 
as  the  number  of  desired  roots  (see  Reference  2).  The  user  should  specify  appropriate  values 
accordingly.  The  shift  point  specified  has  a  significant  effect  on  the  actual  eigensolutions  extracted. 
Accordingly,  depending  upon  the  shift  point  specified,  the  actual  number  of  roots  computed  may  be 
more  or  less  than  the  number  of  desired  roots  specified  in  the  data. 

When  using  the  Inverse  Power  method,  the  EIGR  bulk  data  card  for  each  METHOD  request 
requires  the  specification  of  a  frequency  range,  the  number  of  desired  roots  as  well  as  the  number 
of  estimated  roots  in  the  specified  frequency  range  (see  Reference  2).  The  number  of  estimated  roots 
specified  has  a  significant  effect  on  the  actual  eigensolutions  extracted.  However,  the  user,  in 
general,  will  not  have  an  a  priori  idea  of  the  actual  number  of  eigenvalues  that  may  exist  in  a 
particular  frequency  range.  Accordingly,  the  user  should  use  his  best  judgment  to  specify  this 
number.  A  number  for  the  estimated  roots  that  is  larger  than  the  actual  number  of  roots  in  that  range 
w  ill,  in  general,  yield  a  number  of  lower  frequencies  that  are  outside  the  specified  frequency  range. 

It  should  also  be  noted  here  the  eigensolutions  resulting  from  any  particular  subcase  will  include 
not  only  the  eigensolutions  that  are  computed  in  that  subcase,  but  also  the  eigensolutions  resulting 
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from  all  previous  subcases.  Accordingly,  regardless  of  which  eigenvalue  extraction  method  is 
employed,  the  number  of  desired  roots  specified  on  the  EIGR  card  for  a  particular  subcase  must 
allow  not  only  for  the  roots  that  will  be  computed  by  that  subcase,  but  also  include  the  roots  computed 
by  all  of  the  previous  subcases. 

Output  from  the  Analysis 

As  mentioned  earlier,  the  final  results  from  the  analysis  will  include  the  eigensolutions  for  all 
of  the  subcases  specified  in  the  Case  Control  Deck  setup.  The  DMAP  ALTER  package  given  in 
Appendix  A  also  generates  the  eigenvalue  summary  table  and  the  eigenvalues  for  all  of  the  subcases 
of  the  analysis.  If  the  user  so  desires,  he  can  suppress  any  of  these  intermediate  results  by 
commenting  out  the  OFP  modules  corresponding  to  those  subcases  (see  Appendix  A).  Also,  for 
every  subcase,  the  program  indicates  the  number  of  roots  from  all  previous  subcases  that  are  included 
in  that  subcase. 


EXAMPLES 

Two  examples  were  set  up  to  illustrate  the  procedure  discussed  above.  Details  are  given  below. 
All  of  the  runs  were  made  on  RPK’s  CRAY  version  of  NASTRAN. 

Example  1 

The  standard  NASTRAN  Demonstration  Problem  No.  D 10-02-1 A  (see  Reference  3)  was 
selected  for  this  example.  This  problem  employs  the  Givens  method  for  eigenvalue  extraction.  The 
cyclic  frequencies  obtained  for  this  case  are  presented  in  Table  1. 

This  problem  was  then  modified  to  use  the  Inverse  Power  method  and  eigensolutions  for  two 
frequency  ranges  (500.0  -  1000.0  hertz  and  20000.0  -  30000.0  hertz)  were  requested  using  the 
procedure  described  above.  The  input  data  setup  is  given  in  Appendix  B. 

The  cyclic  frequencies  obtained  for  this  case  are  presented  in  Table  2.  It  can  be  seen  that  these 
frequencies  are  subsets  of  those  in  Table  1. 


Example  2 

A  variation  of  NASTRAN  Demonstration  Problem  No.  D03-08-1A  (without  any  SUPORT 
data)  (see  Reference  3)  was  selected  for  this  example.  This  problem  also  employs  the  Givens  method 
for  eigenvalue  extraction.  The  cyclic  frequencies  obtained  for  this  case  are  presented  in  Tabic  3. 

This  problem  was  then  modified  to  use  the  FEER  method  and  eigensolutions  around  three  shift 
points  or  frequencies  (100.0  hertz,  1500.0  hertz  and  5300.0  hertz)  were  requested  using  the 
procedure  described  above.  The  input  data  setup  is  given  in  Appendix  C. 

The  cyclic  frequencies  obtained  for  this  case  are  presented  in  Table  4.  It  can  be  seen  that  these 
frequencies  are  subsets  of  those  in  Table  3. 
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The  results  of  the  above  examples  clearly  demonstrate  the  validity  and  usefulness  of  the  pro¬ 
posed  method. 


ACKNOWLEDGMENT 

The  procedure  described  above  was  developed  as  part  of  RPK’s  continuing  support  of  its  CRAY 
version  of  NASTRAN  at  NASA’s  Marshall  Space  Flight  Center  in  Huntsville,  Alabama.  The 
authors  are  thankful  to  Mr.  David  Christian  of  NASA-MSFC  for  triggering  the  study  that  made  this 
work  possible. 


CONCLUDING  REMARKS 

A  novel  procedure  for  obtaining  eigenvalues  and  eigenvectors  for  multiple  frequency  ranges  in 
a  single  NASTRAN  execution  is  presented.  The  scheme  is  applicable  to  normal  modes  analyses 
employing  the  FEER  and  Inverse  Power  methods  of  eigenvalue  extraction.  The  procedure  is 
illustrated  by  examples.  The  procedure  should  be  particularly  helpful  in  large  problems  with  widely 
separated  frequency  ranges. 
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APPENDIX  A 


DMAP  ALTERS  for  Obtaining  Eigensolutions  for  Multiple 
Frequency  Ranges  in  a  Single  NASTRAN  Execution 


%  %  5?:  %  %  %  %  &  £  %  #  :*<  Jk  %  %  %  ^  J?:  :k  ^  ;}:  ;fc  ;fc  ^  ;*c  &  5}c  :jc  %  %  sk  #  ❖  &  3j«  ;*<  %  %  ^  Jfc  %  Jfc  ^  &  &  %  Jk  &  ^  ^  &  % 

$  THE  FOLLOWING  ALTERS  ARE  FOR  DISPLACEMENT  RIGID  FORMAT  3 
$  (NORMAL  MODES  ANALYSIS).  SIMILAR  ALTERS  WILL  WORK  FOR 
$  OTHER  RIGID  FORMATS  THAT  INVOLVE  REAL  EIGENVALUE 
$  EXTRACTION. 

$ 

$  NOTE  THAT,  FOR  EVERY  SUBCASE  BEYOND  THE  FIRST  SUBCASE,  THE 
$  ALTERS  BELOW  INVOLVE  A  PAIR  OF  OFP  AND  READ  MODULES. 

$ 

$  INSERT  AFTER  THE  READ  MODULE  IN  THE  RIGID  FORMAT 
$ 


INSERT  READ  $  ON  RPK-SUPPORTED  VERSIONS 

$  USE  ALTER  70  $  ON  1989  COSMIC-SUPPORTED  VERSIONS 

$  USE  OF  ALTER  70  $  IS  ALSO  PERMITTED  ON  RPK-SUPPORTED  VERSIONS 

$ 

$  USE  THE  FOLLOWING  OFP  STATEMENT  TO  REQUEST  THE  EIGENVALUE 
$  SUMMARY  TABLE  AND  THE  EIGENVALUES  THAT  ARE  AUTOMATICALLY 
$  OBTAINED  BY  THE  RIGID  FORMAT  FOR  THE  FIRST  SUBCASE 
$ 


OFP  OEIGS,LAMA„„//S,N,CARDNO  $ 

$ 

$  COMPUTE  THE  EIGENSOLUTIONS  FOR  THE  SECOND  SUBCASE 
$  (THE  LAST  PARAMETER  2  IN  THE  FOLLOWING  READ  MODULE 
$  INDICATES  THAT  THE  ANALYSIS  IS  FOR  THE  SECOND  SUBCASE) 

$ 

READ  KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 
*MODES*/S,N,NEIGV/2  $ 

$ 

$  USE  THE  FOLLOWING  OFP  STATEMENT  TO  REQUEST  THE  EIGENVALUE 
$  SUMMARY  TABLE  AND  THE  EIGENVALUES  FOR  THE  SECOND  SUBCASE 
$  (THE  RESULTS  WILL  INCLUDE  THE  RESULTS  FOR  THE  FIRST  SUBCASE) 
$ 

OFP  OEIGS,LAMA„„//S,N,CARDNO  $ 

$ 

$  COMPUTE  THE  EIGENSOLUTIONS  FOR  THE  THIRD  SUBCASE 
$  (THE  LAST  PARAMETER  3  IN  THE  FOLLOWING  READ  MODULE 
$  INDICATES  THAT  THE  ANALYSIS  IS  FOR  THE  THIRD  SUBCASE) 

$ 
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APPENDIX  A 
(Continued) 


READ  KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 
*MODES  */S  ,N,NEIG  V/3  $ 

$ 

$  USE  THE  FOLLOWING  OFP  STATEMENT  TO  REQUEST  THE  EIGENVALUE 
$  SUMMARY  TABLE  AND  THE  EIGENVALUES  FOR  THE  THIRD  SUBCASE 
$  (THE  RESULTS  WILL  INCLUDE  THE  RESULTS  FOR  THE  FIRST  AND 
$  SECOND  SUBCASES) 

$ 

OFP  OEIGS,LAMA„„//S,N,CARDNO  $ 


$ 

$  COMPUTE  THE  EIGENSOLUTIONS  FOR  THE  LAST  SUBCASE 
$  (THE  LAST  PARAMETER  n  IN  THE  FOLLOWING  READ  MODULE 
$  SHOULD  BE  AN  INTEGER  VALUE  CORRESPONDING  TO  THE  LAST 
$  SUBCASE,  INDICATING  THAT  THE  ANALYSIS  IS  FOR  THE  LAST  SUBCASE) 
$ 

$  THE  FINAL  RESULTS  (WHICH  INCLUDE  THE  RESULTS  FOR  ALL  OF  THE 
$  SUBCASES)  ARE  AUTOMATICALLY  OUTPUT  BY  THE  RIGID  FORMAT 
$ 

READ  KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 
*MODES*/S,N,NEIGV/n  $ 

$ 

CASE  CASECC,/CASEXX/*TRAN*  $ 

EQUIV  CASEXX,CASECC  $ 

ENDALTER  $ 
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APPENDIX  B 


Input  Data  Setup  for  NASTRAN  Demonstration  Problem 
No.  D 10-02-1 A  Modified  to  Use  the  Procedure  Described 

in  the  Paper 


ID  ... 


READFILE  alters 


CEND 


SUBCASE  10 

METHOD  =  100 
SUBCASE  20 

METHOD  =  200 


BEGIN  BULK 


EIGR,  1 00, INV, 500.0, 1 000.0,5, 5,,,+EIG  1 
+EIGLMAX 

EIGR, 200, IN  V, 20000.0, 30000.0, 10, 10,„+EIG2 
+EIG2,MAX 


ENDDATA 
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APPENDIX  C 


Input  Data  Setup  for  NASTRAN  Demonstration  Problem 
No.  D03-08-1A  (Without  Any  SUPORT  Data)  Modified  to 
Use  the  Procedure  Described  in  the  Paper 


ID... 


READFILE  alters 


CEND 


SUBCASE  11 

METHOD  =1001 
SUBCASE  21 

METHOD  =  2001 
SUBCASE  31 

METHOD  =  3001 


BEGIN  BULK 


EIGR,1001,INV,100.0„J5„,+EIG1 

+EIGLMAX 

EIGR, 2001, INV, 1500.0,,, 10,„+EIG2 
+EIG2,MAX 

EIGR, 3001, INV, 5300.0,,, 12,„+EIG3 
+EIG3.MAX 


ENDDATA 
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TABLE  1 


Natural  Frequencies  for  NASTRAN 
Demonstration  No.  D10-02-1A 


Mode  No. 

Cyclic  Frequency  (Hz) 

1 

3.996147E+01 

2 

2.3641 37E+02 

3 

2.504423E+02 

4 

7.01401 1E+02 

5 

7.034198E+02 

6 

1.153105E+03 

7 

1.375429E+03 

8 

1.574398E+03 

9 

1.956923E+03 

10 

2.277239E+03 

11 

2.291 26 1E+03 

12 

2.569 183E+03 

13 

2.783842E+03 

14 

2.929954E+03 

15 

3.00392 1E+03 

16 

3.41 1562E+03 

17 

4.786588E+03 

18 

6.412708E+03 

19 

8.291552E+03 

20 

1.030759E+04 

21 

1.371895E+04 

22 

1.657422E+04 

23 

2.0091 81E+04 

24 

2.424947  E+04 

25 

2.913932E+04 

26 

3.485085E+04 

27 

4.136255E+04 

28 

4.82888 1E+04 

29 

5.437856E+04 

30 

6.8054 13E+04 
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TABLE  2 


Natural  Frequencies  for  NASTRAN  Demonstration 
No.  D10-02-1A  Modified  to  Use  the  Procedure 
Described  in  the  Paper 


Computed 

Mode 

No. 

Corresponding 
Mode  No. 
in  Table  1 

Cyclic 
Frequency 
(Hz)  ' 

Subcase  in  Which 
Computed 
(See  Appendix  B) 

1 

2 

2.364 137E+02 

First 

2 

3 

2.504423E+02 

First 

3 

4 

7.01401 1E+02 

First 

4 

5 

7.034 198E+02 

First 

5 

22 

1.657422E+04 

Second 

6 

23 

2.009181E+04 

Second 

7 

24 

2.424947E+04 

Second 

8 

25 

2.913932E+04 

Second 

9 

26 

3.485085E+04 

Second 
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TABLE  3 


Natural  Frequencies  for  NASTRAN 
Demonstration  No.  D03-08-1A 
(Without  Any  SUPORT  Data) 


ode  No. 

Cyclic  Frequency  (Hz) 

1 

7.341481E-04 

2 

4.168523E-04 

3 

6.790643E-05 

4 

1.505714E-04 

5 

2.765971E-04 

6 

6.434454E-04 

7 

2.987344E+00 

8 

3.372945E+00 

9 

2.447569E+01 

10 

2.682217E+01 

11 

6.154903E+01 

12 

7.034309E+01 

13 

1.133579E+02 

14 

1.1 7453 1E+02 

15 

1.646037E+02 

16 

2.902883E+02 

17 

2.905903E+02 

18 

4.508246E+02 

19 

4.515298E+02 

20 

6.945689E+02 

21 

6.953026E+02 

22 

8.685333E+02 

23 

9.752645E+02 

24 

9.752699E+02 

25 

1.343024E+03 

26 

1.344180E+03 

27 

1.466313E+03 

28 

1.569195E+03 

29 

1.912914E+03 

30 

1.918371E+03 

31 

2.02598 1E+03 

32 

2.446537E+03 

33 

2.446840E+03 

34 

2.458828E+03 

35 

2.742675E+03 

36 

2.971 822E+03 

37 

3.918531E+03 

38 

3.918534E+03 

39 

4.451 85 1E+03 

40 

6.261906E+03 

41 

8.528055E+03 

TABLE  4 


Natural  Frequencies  for  NASTRAN  Demonstration 
No.  D03-08-1A  (Without  Any  SUPORT  DATA)  Modified 
to  Use  the  Procedure  Described  in  the  Paper 


Computed 

Mode 

No. 

Corresponding 
Mode  No. 
in  Table  3 

Cyclic 
Frequency 
(Hz)  ' 

Subcase  in  Which 
Computed 
(See  Appendix  C) 

1 

9 

2.447569E+01 

First 

2 

10 

2.682217E+01 

First 

3 

11 

6.154903E+01 

First 

4 

12 

7.034309E+01 

First 

5 

13 

1.133579E+02 

First 

6 

14 

1.174531E+02 

First 

7 

25 

1.343024E+03 

Second 

8 

26 

1.344180E+03 

Second 

KB 

27 

1.4663 13E+03 

Second 

28 

1.569195E+03 

Second 

39 

4.451 85 1E+03 

Third 

40 

6.261906E+03 

Third 

RANDOM  VIBRATION  ANALYSIS 


OF  SPACE  FLIGHT  HARDWARE  USING  NASTRAN 

S.  K.  Thampi  and  S.  N.  Vidyasagar 
GE  Government  Services 
'  1 050  Bay  Area  Blvd. 

Houston,  TX  77058 


During  liftoff  and  ascent  flight  phases,  the  Space  Transportation  System  (STS) 
and  payloads  are  exposed  to  the  random  acoustic  environment  produced  by 
engine  exhaust  plumes  and  aerodynamic  disturbances.  The  analysis  of 
payloads  for  randomly  fluctuating  loads  is  usually  carried  out  using  the  Miles’ 
relationship.  This  approximation  technique  computes  an  equivalent  load  factor 
as  a  function  of  the  natural  frequency  of  the  structure,  the  power  spectral  density 
of  the  excitation,  and  the  magnification  factor  at  resonance.  Due  to  the 
assumptions  inherent  in  Miles’  equation,  random  load  factors  are  often  over¬ 
estimated  by  this  approach.  In  such  cases,  the  estimates  can  be  refined  using 
alternate  techniques  such  as  time  domain  simulations  or  frequency  domain 
spectral  analysis.  This  paper  describes  the  use  of  NASTRAN  to  compute  more 
realistic  random  load  factors  through  spectral  analysis.  The  procedure  is 
illustrated  using  Spacelab  Life  Sciences  (SLS-1 )  payloads  and  certain  unique 
features  of  this  problem  are  described.  The  solutions  are  compared  with  Miles' 
results  in  order  to  establish  trends  at  over  or  under  prediction. 
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INTRODUCTION 


During  the  past  decade,  the  U.S.  Spacelab  program  has  made  significant 
contributions  to  the  advancement  of  space  exploration  and  research.  The 
Spacelab  is  a  reusable  laboratory  that  is  carried  in  the  cargo  bay  of  the  Space 
Shuttle.  Experiments  in  several  different  disciplines  such  as  astronomy,  life 
sciences,  and  material  science  are  accommodated  in  this  modular  laboratory 
for  various  shuttle  missions.  The  module  also  contains  utilities,  computers,  and 
work  areas  to  support  the  experiments.  The  experiment  hardware  is  mounted  in 
instrument  racks  located  on  either  side  of  the  module,  in  overhead  lockers,  or  in 
the  center  aisle,  as  shown  in  Figure  1. 

During  liftoff  and  ascent  flight  events,  the  Shuttle  and  its  payload  are  exposed  to 
the  acoustic  environment  produced  by  engine  exhaust  plumes  and 
aerodynamic  disturbances.  Random  vibrations  are  created  by  the  response  of 
the  module  shell  to  the  acoustic  noise  inside  the  cargo  bay.  The  vibrations  of 
the  shell  are  transmitted  through  the  support  structures  (racks,  mounting  frames, 
etc.)  to  the  payload  equipment.  The  vibration  levels  that  the  equipment  has  to 
withstand  depend  on  its  own  dynamic  characteristics  and  its  location  inside  the 
Spacelab.  The  equipment  and  its  structural  interfaces  must  be  analyzed  for 
these  random  loads  in  order  to  ensure  the  integrity  and  flight  worthiness  of  the 
system. 

The  analysis  of  flight  hardware  for  random  loads  often  relies  on  approximate 
formulations  like  the  Miles'  relation  (ref.  1)  to  generate  limit  load  factors  for 
structural  design.  Due  to  the  assumptions  inherent  in  Miles’  equation,  the 
random  vibration  criteria  developed  through  this  approach  tend  to  be  over¬ 
conservative.  In  such  cases,  the  results  can  be  refined  using  alternate  analysis 
techniques  like  time  domain  simulations  or  frequency  domain  spectral  analysis. 
This  paper  describes  the  use  of  NASTRAN  to  perform  spectral  analysis  to 
establish  more  realistic  design  loads.  The  procedure  is  illustrated  using  the 
Neck  Chamber  Pressure  System  (NCPS)  assembly  which  will  be  flown  on  the 
Spacelab  Life  Sciences  (SLS-1)  mission. 


ANALYSIS  BASED  ON  MILES'  EQUATION 


For  a  lightly  damped  single-degree-of-freedom  (SDOF)  oscillator  subjected  to 
random  excitation  through  its  base  motion,  Miles'  relation  is  used  as  follows. 

N,  =  3  -\jj  f„  A  Q  (1) 

Nr  is  the  limit  random  load  factor  (g  units),  fn  is  the  resonance  frequency  (Hz)  of 
the  SDOF  system,  A  is  the  power  spectral  density  (g2/Hz)  of  base  acceleration  at 
the  resonance  frequency,  and  Q  is  the  dynamic  magnification  factor  at 
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resonance.  The  quantity  under  the  square  root  represents  the  mean  square 
acceleration  response  and  the  factor  3  indicates  the  3-sigma  probability  of 
occurrence,  i.e.,  the  probability  of  exceeding  Nr  is  0.26%. 

The  mean  square  response  is  the  area  under  the  response  spectral  density 
curve  and  is  given  by 


<u2>  =  J  S(co)  |H(co)|2  dco  (2) 

where  H(co)  is  the  transfer  function  of  the  system,  S  (co)  is  the  base  acceleration 

spectral  density  function,  and  co  is  the  frequency  in  radians.  The  derivation  of 
Miles'  relation  is  based  on  the  following  simplifying  assumptions  for  evaluating 
the  integral  in  Eqn  (2). 

1 )  The  actual  spectral  density  of  base  excitation,  S  (co),  is  a  slowly  varying 
function  in  the  vicinity  of  resonance.  It  can  be  conveniently  approximated  by  a 
constant  or  white-noise  spectral  density,  A  =  S  (con),  for  computing  the  mean 
square  response. 

2)  Only  the  excitation  energy  contained  within  the  system’s  bandwidth  is 
transmitted.  The  rest  is  filtered  away  by  the  system  and  does  not  contribute  to 
the  mean  square  response. 

These  assumptions  are  valid  in  the  case  of  lightly  damped  systems  with 
damping  factor  £,«}.  For  such  systems,  the  function  |H(co)|2  is  very  sharply 

peaked  at  co  =  con ,  and  reduces  to  half  its  peak  value  at  a  short  distance,  2^con, 
on  either  side  of  the  peak.  This  distance,  called  the  half  power  bandwidth,  is 
very  narrow  for  lightly  damped  systems.  With  the  assumptions  mentioned 
above,  the  integral  in  Eqn.  2  can  be  approximated  by  a  rectangular  area  with 
base  equal  to  the  bandwidth  and  height  equal  to  the  product  of  the  constant 
value  of  excitation  spectral  density  and  the  peak  value  of  transmittancy.  This 
gives 


<u2>  ~  I?  fnAQ  (3) 

from  which  Miles'  relation  follows. 

In  order  to  use  Miles'  relation  for  the  analysis  of  flight  hardware,  the  natural 
frequency  of  the  equipment  is  first  determined  through  analysis  or  test.  The 
input  random  excitation  spectrum  for  the  equipment  is  then  determined  as  a 
function  of  its  location,  its  mounting  configuration,  and  its  mass.  The  input 
excitation  spectrum  has  been  established  using  data  from  previous  flight  or 
ground  tests  and  is  provided  in  the  Spacelab  Payload  Accomodation  Handbook 
(ref.  2).  The  spectral  density  at  the  resonance  frequency  of  the  component  is 

found  from  this  data.  The  dynamic  magnification  factor  at  resonance  Q  =  1/2  £, 
is  indicative  of  the  system  damping  and  is  determined  experimentally.  For 
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components  mounted  on  isolators,  Q  is  determined  from  the  manufacturer’s 
data  on  the  isolator  mounting.  These  values  are  substituted  in  Eqn  (1)  to  obtain 
the  design  random  load  factor,  Nr.  As  the  random  vibration  environment  occurs 
simultaneously  with  other  load  environments  during  various  mission  phases, 
the  estimated  values  of  Nr  are  combined  with  the  appropriate  quasi-static, 
thermal,  pressure,  and  crew-induced  load  factors  to  generate  design  load  cases 
for  component  analysis. 

Due  to  the  assumptions  inherent  in  Miles'  relation  and  the  idealization  of  the 
component  as  a  single-degree-of-freedom  resonator,  the  computed  random 
load  factors  will  be  approximate.  They  tend  to  be  overly  conservative, 
especially  when  the  natural  frequency  of  the  system  is  close  to  the  peak 
frequency  of  the  excitation  spectrum.  When  the  predicted  random  loads  are 
unreasonably  high,  they  lead  to  difficult  design  problems  and  alternate 
approaches  are  necessary  to  refine  the  random  load  estimates. 


ANALYSIS  BASED  ON  SPECTRAL  ANALYSIS 


The  dynamic  behaviour  of  large  structural/mechanical  systems  can  be 
adequately  predicted  only  by  multi-degree-bf-freedom  (MDOF)  models.  For 
linear  MDOF  systems,  the  dynamic  characteristics  are  specified  by  a  matrix  of 
transfer  functions,  [H],  whose  elements  Hjk  represent  the  ratio  of  steady-state 
response  at  point  j  to  a  sinusoidal  excitation  at  point  k.  For  displacement 
response 


[H(co)J  =  [-[M]  co2  +  i  [B]  co  +  [K]]  -i  (4) 

where  [M],  [B],  and  [K]  represent  the  mass,  damping,  and  stiffness  matrices  of 
the  discrete  model  and  co  is  the  excitation  frequency. 

The  response  of  linear  MDOF  systems  subjected  to  random  excitation  can  be 
computed  using  well-established  spectral  analysis  techniques.  According  to 
the  theory  of  random  vibrations,  the  response  of  a  linear  system  with  transfer 
function  [H(co)j,  subjected  to  a  stationary  random  load  {P(t)}>  is  given  by 

[Suu(co)j  =  [H(co)j  [Spp(co)]  [H*(co)]T  (5) 

where  [Suu  (co)]  and  [Spp(co)]  are  the  matrices  of  response  and  excitation  spectral 
density  functions  and  *  and  T  represent  the  complex  conjugate  and  transpose 
operations,  respectively.  These  matrices  will  have  real  auto-spectral  density 
functions  as  their  diagonal  elements  and  complex  cross-spectral  density 
functions  as  their  off-diagonal  elements.  By  integrating  the  area  under  the 
response  spectral  density  curve,  the  mean  square  response  at  any  nodal  point 
in  the  model  can  be  obtained.  The  foregoing  development  for  mean-square 
displacement  response  can  be  generalized  to  provide  mean-square  values  for 
other  response  quantities  such  as  velocity,  acceleration,  or  stress.  It  is  only 
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necessary  to  replace  the  transfer  function  matrix  for  displacement  by  the 
corresponding  transfer  function  matrix  for  the  desired  response.  The 
generalization  also  applies  to  the  excitation  which  may  be  a  point  force,  a 
loading  condition  (i.e.,  an  ensemble  of  applied  forces  that  are  perfectly 
correlated),  or  enforced  motion. 

The  analysis  of  flight  hardware  subjected  to  random  excitation  can  be  carried 
out  using  the  spectral  analysis  features  of  NASTRAN.  A  finite  element  model  of 
the  component  is  created  which  can  accurately  represent  all  its  dominant 
modes  in  the  excitation  frequency  range.  The  response  calculations  are  carried 
out  in  two  separate  functional  modules.  First,  the  transfer  function  of  the  system 
corresponding  to  the  desired  response  is  computed  in  the  Frequency  Response 
module  and  then,  the  power  spectral  densities  and  other  response  statistics  are 
computed  in  the  Random  Analysis  module.  The  direct  or  modal  superposition 
approaches  can  be  used  to  perform  the  frequency  response  analysis.  For  each 
excitation  source,  pk,  the  nodal  response,  Uj,  is  determined  at  a  series  of  user 

specified  frequencies,  coj.  The  ratio  of  output  to  input  represents  the  transfer 

function  element,  Hjk(coj ).  This  is  determined  for  each  excitation  source  and  the 
transfer  function  matrix,  [H],  is  assembled  from  the  results. 

In  NASTRAN,  random  vibration  analysis  is  treated  as  a  data  reduction 
procedure  that  is  applied  to  the  results  of  frequency  response  analysis.  The 
inputs  to  the  random  analysis  module  are  the  frequency  responses  of  desired 
output  quantities  due  to  different  load  sources  and  the  auto-  and  cross-spectral 
densities  of  these  random  load  sources.  Each  load  source  is  referred  to  by  a 
separate  subcase  in  the  case  control  deck  and  their  spectral  densities  are 
specified  as  tabular  functions  of  frequency  in  bulk  data  cards.  If  the  sources  are 
statistically  uncorrelated,  only  the  auto-spectral  densities  need  be  defined.  The 
power  spectral  densities  of  response  are  calculated  using  Eqn  (5)  and  the  root- 
mean-square  (rms)  response  is  evaluated  by  numerically  integrating  the  area 
under  the  spectral  density  curve.  The  results  are  printed  and  plotted  for 
specified  degrees  of  freedom  of  the  model. 

As  mentioned  earlier,  the  random  excitation  applied  to  the  structure  could  be  a 
force,  an  enforced  motion,  or  some  other  general  form  of  excitation.  In  the  case 
of  Spacelab  payloads,  the  random  excitation  is  specified  in  terms  of  an 
acceleration  spectrum  applied  at  the  structural  support  points.  The  "large  mass" 
approach  may  be  used  to  simulate  this  loading  condition.  This  involves 
lumping  a  fictitious  large  mass,  Ma,  at  the  degree  of  freedom  in  which  the 
acceleration  is  to  be  enforced.  An  applied  force  equal  to  Ma  times  the  required 
acceleration  is  also  prescribed  for  that  degree  of  freedom.  The  inertia  force  is 
made  so  dominant  through  this  operation  that  the  resulting  acceleration  is  very 
close  to  the  required  value. 
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ILLUSTRATIVE  EXAMPLE 


The  development  of  random  vibration  load  factors  for  flight  hardware  is 
illustrated  in  this  section  using  the  Neck  Chamber  Pressure  System  assembly 
(NCPS).  This  LM  Sciences  Laboratory  equipment  item  is  used  to  study  the 
effect  of  weightlessness  on  human  cardiovascular  control  mechanisms  arid  will 
be  flown  on  a  future  Space  Shuttle  mission.  The  NCPS  assembly  is  composed 
of  an  experiment  enclosure  which  houses  several  components  including  a 
central  processing  unit,  a  motor  control  unit,  two  motor  driven  bellows,  and  a 
pressure  gauge  (Fig.  2).  The  assembly  is  mounted  in  an  experiment  rack  using 
two  support  rails  which  are  attached  to  the  front  and  rear  rack  posts  on  either 
side.  The  front  panel  of  the  enclosure  is  bolted  to  the  front  rack  post  flanges  at 
eight  locations.  The  whole  assembly  weighs  48  pounds,  and  the  installation-kit 
including  the  slides,  fittings,  and  fasteners  weighs  an  additional  five  pounds. 

A  finite  element  model  of  the  NCPS  assembly  is  constructed  using  mostly  plate 
(CQUAD2  and  CTRIA2)  and  bar  (CBAR)  elements.  The  model  has  a  total  of  206 
grid  points  and  192  structural  elements.The  masses  of  internal  components  are 
lumped  at  the  respective  centers  of  gravity,  and  stiff  bar  elements  are  used  to 
connect  them  to  the  attachment  points.  The  fasteners  are  modelled  using  rigid 
elements.  Eigenanalysis  was  performed  on  the  model' with  free  boundary 
conditions  to  verify  that  the  model  has  six  rigid  body  modes.  The  analysis  is 
repeated,  with  the  rack-to-component  interface  points  appropriately 
constrained,  in  order  to  determine  the  flexible  modes  of  the  component.  The, 
first  twenty  frequencies  of  the  constrained  model  are  shown  in  Table  1 .  An 
inspection  of  the  modal  deformation  plots  and  mass  participation  factors  shows 
that  the  first  system  mode  in  X  direction  is  80  Hz.  The  power  spectral  density  of 
the  input  excitation  corresponding  to  this  frequency  is  found  to  be  0.02  g2/Hz 
from  Figure  3.  For  a  conservative  estimate  of  Q  =  10  for  the  dynamic 
magnification  factor,  Miles'  relation  (Eqn  1)  yields 

Nrx  =  15.04  g  units  (6) 

Random  load  factors,  N,y  and  Nr2)  are  computed  in  a  similar  manner. 

Random  load  factors  can  also  be  determined  through  spectral  analysis.  The 
computation  of  Nrx  is  described  here  for  comparison  with  the  Miles'  approach. 
The  transfer  function  of  the  system  is  first  determined  by  applying  a  unit 
sinusoidal  load  in  the  X  direction  at  each  point  where  the  NCPS  interfaces  with 
the  rack  structure.  The  load  is  applied  through  the  DLOAD  and  RLOAD  cards 
which  in  turn  refers  to  DELAY,  DPHASE,  DAREA,  and  TABLED  cards.  For  a 
constant  phase,  unit  sinusoidal  input,  the  DELAY  and  DPHASE  cards  may  be 
omitted  and  a  unit  value  specified  for  all  frequencies  on  the  TABLED  card.  The 
input  acceleration  spectrum  (Fig.  3)  is  specified  through  RANDPS  and  TABRND 
cards.  This  random  acceleration  is  enforced  at  all  the  interface  points  in  the  X 
direction  using  the  large  mass  approach.  A  fictitious  mass,  Ma,  about  1000 
times  larger  than  the  existing  grid  point  mass,  is  lumped  at  the  interface  X 
degree  of  freedom  using  a  CMASS2  card  and  an  equal  value  is  specified  in  the 
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corresponding  DAREA  load  card.  This  produces  the  desired  acceleration  which 
can  be  verified  by  plotting  the  acceleration  spectrum  at  the  input  points  and 
comparing  it  with  Figure  3.  If  the  agreement  is  not  adequate,  the  value  of  Ma  is 
increased  until  a  good  match  is  obtained. 

The  direct  and  modal  solution  techniques  are  used  to  carry  out  the  response 
calculations.  The  relative  efficiency  of  the  two  approaches  depends  on  several 
factors  including  the  number  of  modes  retained  in  the  analysis  and  the  number 
of  response  frequencies.  The  selection  of  these  parameters  always  represents 
a  compromise  between  accuracy  and  efficiency.  The  frequencies  chosen  for 
response  computations  should  have  good  resolution  in  the  vicinity  of  system 
resonances  in  order  to  obtain  reliable  estimates  of  rms  response.  A  total  of  150 
frequencies  in  the  0  to  500  Hz  range  are  used  for  this  analysis.  The 
specification  of  damping  properties  in  the  direct  and  modal  formulations  are 
somewhat  different.  In  the  direct  formulation,  structural  damping  proportional  to 
the  stiffness  matrix  terms  is  specified  both  on  the  material  data  cards  and  as  an 
overall  uniform  damping  factor  on  the  PARAM  G  data  card.  In  the  modal 
formulation,  the  damping  factor  is  specified  as  a  tabular  function  of  frequency 
through  the  SDAMPING  and  TABDMP1  cards  to  represent  the  variation  of 
structural  damping  for  different  modes.  A  damping  factor  of  0.1  is  used  for  this 
analysis  in  order  to  be  consistent  with  Q  =  10  used  in  Miles'  relation. 

The  set  of  output  points  at  which  the  response  power  spectrum  and  rms  values 
are  to  be  recovered  must  be  chosen  judiciously.  The  selection  can  be  based  on 
the  same  criteria  used  for  choosing  an  ASET  (analysis  set  of  dynamic  degrees 
of  freedom);  namely,  the  points  should  be  uniformly  dispersed  throughout  the 
structure  and  should  include  all  large  mass  items.  A  set  of  12  output  points 
were  chosen  for  the  NCPS  and  the  rms  acceleration  response  at  these  points  is 
computed  (Table  2).  When  these  values  are  averaged  and  the  3-sigma 
probability  of  occurrence  criteria  is  applied,  one  obtains 

No?  =  10.67  g  units  (7) 

This  is  almost  30%  less  than  the  value  predicted  using  Miles'  relation. 

The  rms  response  calculated  using  the  direct  and  modal  solution  techniques  is 
summarized  in  Table  2.  For  the  same  number  of  response  frequencies,  the 
modal  solution  required  350  seconds  of  CPU  time  with  20  modes  included,  450 
seconds  with  40  modes,  600  seconds  with  60  modes,  whereas  the  direct 
solution  required  13160  seconds.  The  response  spectrum  at  selected  output 
points  on  the  model  is  shown  in  figures  4,  5,  and  6.  The  solution  obtained  with 
20  modes  is  clearly  inadequate  while  the  plots  obtained  using  60  modes  and 
the  direct  approach  are  virtually  indistinguishable.  The  spectra 
characteristically  peak  at  80  Hz  which  corresponds  to  the  first  X  mode  of 
vibration.  The  effects  of  higher  modes  can  also  be  seen  in  the  spectral  plots. 
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An  alternate  method  of  estimating  random  vibration  load  factors  for  design  and 
analysis  of  Spacelab  payloads  is  presented.  This  method,  based  on  spectral 
analysis,  yields  more  refined  random  load  estimates  at  the  expense  of  being 
computationally  more  intensive  than  the  Miles'  approach.  The  computational 
effort  can  be  reduced  by  using  the  modal  formulation  rather  than  the  direct 
formulation  for  analysis.Significant  reductions  can  be  obtained  in  the  random 
load  estimates  using  this  method.  While  the  actual  reduction  depends  upon  the 
payload  configuration  being  analyzed,  reductions  of  20  to  30%  are  typical.  This 
method  could  be  used  to  resolve  difficult  design  problems  owing  to 
unreasonably  high  random  load  predictions  by  the  Miles'  relation. 


REFERENCES 


1.  Harris,  C.  M.  and  Crede,  C.  E.,  "Shock  and  Vibration  Handbook,"  2nd 
Edition,  McGraw  Hiii  Book  Company,  Inc.,  New  York,  1976,  p‘p.24-28. 

2.  Spaceiab  Payload  Accomodation  Handbook:  Appendix  B,  SLP/2104-2, 
NASA  Marshall  Space  Flight  Center,  October  1989. 


98 


TABLE  I 


EIGENVALUE  ANALYSIS  SUMMARY 


MODE 

NO. 

EIGENVALUE 

RADIAN 

FREQUENCY 

CYCLIC 

FREQUENCY 

1 

5 . 606682E+04 

2 . 367843E+02 

3 . 768539E+01 

2 

8.630755E+04 

2 . 937815E+02 

4 . 675677E+01 

3 

9 . 900016E+04 

3 . 146429E+02 

5 . 007697E+01 

4 

1. 228704E+05 

3 . 505287E+02 

5 . 578837E+01 

5 

1 . 483150E+05 

3 . 851169E+02 

6 . 129326E+01 

6 

1. 896878E+05 

4 . 355316E+02 

6 . 931701E+01 

7 

2 . 301971E+05 

4 . 797886E+02 

7 . 636072E+01 

8 

2 . 535269E+05 

5 . 035146E+02 

8 . 013683E+01 

9 

3 . 989269E+05 

6 . 316066E+02 

1 . 005233E+02 

10 

4 . 230480E+05 

6 . 504214E+02 

1 . 035178E+02 

11 

4 . 433938E+05 

6 . 658782E+02 

1.059778E+02 

12 

4 . 683294E+05 

6 . 846459E+02 

1 . 08917  0E+02 

13 

5 . 100052E+05 

7 . 141465E+02 

1 . 136599E+02 

14 

5 . 928444E+05 

7 . 699639E+02 

1 . 225435E+02 

15 

7 . 732380E+05 

8 . 793395E+02 

1 . 399512E+02 

16 

7 . 972504E+05 

8 . 928888E+02 

1 . 421077E+02 

17 

8 . 693909E+05 

9 . 324113E+02 

1.483979E+02 

18 

9 . 891926E+05 

9 . 945816E+02 

1 . 582926E+02 

19 

1 . 104014E+06 

1 . 050721E+03 

1 . 672274E+02 

20 

1 . 231224E+06 

1 . 109605E+03 

1 . 7 65992E+02 

21 

1 . 492994E+06 

1 . 22 1881E+03 

1 . 944685E+02 
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TABLE  II  RMS  ACCELERATION  RESPONSE  SUMMARY 


GRID 

POINT 

LOCATION 

MODAL  SOLUTION' 

DIRECT 

SOLUTION 

(g) 

- 

20  MODES 

(g) 

40  MODES 

(g) 

60  MODES 

(g) 

5059 

Rear  Left 
Side  Panel 

2.07 

3.50 

3.52 

3.52 

5080 

Rear  Right 
Side  Panel 

2.06 

3.97 

3.99 

3.98 

5056 

Front  Left 
Side  Panel 

2.47 

2.47 

2.49 

2.50 

5077 

Front  Right 
Side  Panel 

2.47 

2.47 

2.49 

2.50 

5076 

Top  Left 
Rear  Panel 

2.08 

4.53 

4.58 

4.58 

5073 

Top  Left 
Front  Panel 

2.24 

2.48 

2.87 

2.88 

5094 

Top  Right 
Front  Panel 

2.24 

2.48 

2.87 

2.88 

5097 

Top  Right 
Rear  Panel 

2.08 

4.89 

4.91 

4.91 

5193 

M/C  Unit 

2.07 

4.33 

4.46 

4.46 

5204 

CPU  Unit 

2.81 

3.89 

3.95 

3.95 

5205 

M/C  Mount 

2.12 

4.03 

4.04 

4.03 

5042 

Front  Panel 

2.47 

2.47 

2.49 

2.50 

Average 

2.26 

3.46 

3.55 

3.56 

100 


TYPICAL  SPACELAB  CONFIGURATION 


SPECTRAL  DENSITY  (  g  /Hz 
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FIG.  4.  ACCELERATION  RESPONSE  SPECTRUM  AT  GRID  POINT  5073 


(  *n/2*>  ) 


J1XSH3Q  TYai03dS 


105 


fHEOOENCY 


OBTAINING  AN  EQUIVALENT  BEAM 


THOMAS  G.  BUTLER 
BUTLER  ANALYSES 


In  modeling  a  complex  structure  I  was  faced  with  a 
component  that  would  have  logical  appeal  if  it  were  modeled  as  a 
beam.  It  was  a  mast  of  a  robot  controlled  gantry  crane.  The 
structure  up  to  this  poi.it  already  had  a  large  number  of  degrees 
of  freedom,  so  the  idea  of  conserving  grid  points  by  modeling  the 
mast  as  a  beam  was  attractive.  I  decided  to  make  a  separate 
problem  of  the  mast  and  model  it  in  three  dimensions  with  plates 
then  extract  the  equivalent  beam  properties  by  setting  up  the 
loading  to  simulate  beam  like  deformations  and  constraints.  The 
results  could  then  be  used  to  represent  the  mast  as  a  beam  in  the 
full  model.  This  seemed  to  be  a  straight  forward  approach,  but 
it  was  sufficiently  challenging  that  it  merited  publishing  a 
paper  on  this  topic. 

The  endeavor  is  to  obtain  the  area  A,  the  area  moments  of 
inertia  II  and  12,  and  torsional  area  moment  of  inertia  J  of  a 
prismatic  beam  that  would  be  an  equivalent  of  the  crane  mast  over 
its  full  length.  The  detailed  model  involved  about  4500  uncon¬ 
strained  degrees  of  freedom.  The  mast  structure  was  essentially 
a  hollow  steel  tube  of  square  section  with  a  cylindrical  indenta¬ 
tion  along  its  length  on  one  surface  only.  Complications  that 
made  it  difficult  to  estimate  equivalent  properties  analytically 
were  the  placement  of  two  types  of  interior  partial  shear  stiff¬ 
eners  at  regular  intervals  along  its  length.  These  two  different 
types  of  shear  stiffeners  alternated  on  opposite  sides  from  each 
other  most  of  the  length.  Thi3  posed  no  difficulty  to  model 
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elastically  in  a  three  dimensional  model.  The  interesting  phase 
is  the  loading  of  the  3-D  model  in  order  to  simulate  beam  action. 

To  put  che  problem  in  perspective,  review  for  che  moraenc , 
the  definition  of  beam  stiffness. 

DEFINITION:  Beam  stiffness  is  the  array  of  forces  pro¬ 
duced  at  the  six  degrees  of  freedom  on  both  ends  when  a 
single  degree  of  freedom  at  one  end  is  deformed  a  unit 
amount  while  enforcing  all  other  eleven  degrees  of 
freedom  at  both  ends  to  be  zero. 

But  the  Bernoulli  Euler  formulation  of  the  beam  as  used  in  finite 
element  analysis  programs  does  not  faithfully  follow  this 
prescription  of  stiffness  to  the  letter.'  For  example,  when  one 
end  is  displaced  a  unit  transversely,  action  is  assumed  to  occur 
in  -  plane  only.  Diagrammatically  the  boundary  conditions  of  the 
centroid  of  the  B.E.  formulation  are  indicated  in  the  sketch. 


Note  that  the  length  remains  invariant,  because  its  transversely 
deformed  end  is  not  constrained  in  the  axial  direction.  In 
effect,  with  this  B.E.  approach,  the  end  position  contracts  when 
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bending  deformation  occurs.  This  is  shown  in  exaggerated  fashion 
in  the  sketch. 


If  the  length  is  not  allowed  to  delorm,  Poisson  deformation  does 
not  occur  and  therefore  needs  no  constraining  force  to  inhibit 
Poisson  deformation.  But  if  the  true  definition  of  beam  stiff¬ 
ness  were  adhered  to  in  the  finite  element  beam,  the  axial  posi¬ 
tions  of  the  ends  would  be  held  to  2;erq  displacement  and  the  beam 
would  lengthen  as  transverse  deformation  occurs.  Such  axial 
stretching  would  result  in  Poisson  contraction  in  both  transverse 
directions.  But  if  transverse  translational  deformations  were 
held  to  zero,  as  the  definition  of  stiffness  demands,  such  con¬ 
straints  would  exert  forces  to  prevent  Poisson  contraction.  For 
instance,  the  transverse  forces  at  the  end  of  a  solid  beam  of 
square  section  with  a  full  set  of  constraints  applied  would 
appear  as  sketched. 
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The  dilemma  now  is  to  try  to  define  what  kind  of  equivalence 
should  be  sought.  If  A,  II,  12,  and  J  were  obtained  with  true 
stiffness  constraints,  would  it  be  proper  to  operate  as  an 
equivalent  beam  according  to  those  entries  on  the  property  card, 
so  as  to  exclude  bending/axial  coupling  even  though  such  action 
was  present  during  the  sample  run?  Or  would  it  be  more  proper  to 
use  only  B.E.  conditions  to  get  the  properties  that  will  used  as 
a  B.E.  beam?  If  the  latter  were  chosen,  the  question  arises  as 
to  how  faithfully  we  would  be  representing  equivalence  to  the 
true  structure.  Having  some  doubts  as  to  how  to  proceed,  I 
modeled  the  constraints  in  two  different  ways;  with  full  end 
constraints  and  with  B.E.  end  constraints  and  compared  the 
results.  The  sketch  shows  the  constraints  imposed  for  the  two 
models.  One  of  the  things  to  consider  in  the  B.E.  simulation  is 
that  the  theory  requires  planes  to  remain  plane  in  bending. 
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The  next  question  is:  After  constraint  forces  are 
measured,  will  it  be  acceptable  to  derive  sectional  properties  by- 
substituting  into  the  formulation  based  strictly  on  B.E.?  That 
is  to  say,  should  the  stiffness  forces  obtained  on  the  left  be 
equated  to  the  B.E.  formulas  on  the  right?  Just  enough  of  the 
matrix  on  the  right  is  shown  to  illustrate  the  problem. 


K 


11 


K 


22 


K17 

K26  K28 


K 


2  C 


K33 

K35 

K39 

K3B 

K44 

K4A 

K53 

K55 

K59 

K5B 

I< 


62 


K66  K68 


CC 


EA/L 


12EI„/L' 

4 


12EI  /L‘ 

y 


GJ/L 


-6EI  /L' 


6EIZ/L 


Not  having  any  reference  to  use  for  the  fully  coupled 
beam  I  chose  to  use  B.E.  formulation  to  evaluate  sectional 
properties  for  both  types  of  modeling. 

The  next  question  is:  After  accepting  B.E.  formulation, 
what  basis  should  be  used  to  reconcile  differences  in  results  of 
the  methods?  The  reconciliation  method  is  to  use  an  estimation 
of  the  computed  value  of  the  section  without  the  shear  panels 
present  as  per  the  dimensions  in  the  sketch. 
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COMPARISON  OF  PROPERTIES  DERIVED  FROM  MODELS  OF  DIFFERENT 
CONSTRAINTS  VS  MANUAL  CALCULATIONS 


SOURCE 

1 

A  | 

11  1 

12  | 

J 

FULL  CONSTRAINTS 

1 

1 

36.72  | 

I 

2,444.47  | 

2,605.14  | 

INVALID 

B.E.  CONSTRAINTS 

1 

1 

1 

1 

36.66  | 

| 

1 

728.63  | 

1 

853.62  | 

600.19 

MANUAL 

1 

1 

1 

35.96  | 

1 

2,541.82  | 

1 

3,517.62  | 

4,710.60 

This  exercise  had  some  unexpected  results.  The  whole 
purpose  of  the  exercise  was  to  get  an  equivalent  beam  by  using  a 
full  3-D  model  instead  of  making  an  analytical  estimate  because 
of  the  uncertainty  in  being  able  to  represent  the  effect  of  the 
partial  shear  panels  correctly.  One  expects  that  the  effect  of 
the  shear  panels  is  to  stiffen  the  steel  tube,  but  the  3-D  re¬ 
sults  showed  less  stiffness  than  the  manual  check  which  neglected 
the  panels.  Why? 

In  going  back  to  examine  the  axial  displacements  in  the 
3-D  model  using  the  B.E.  constraints,  it  indicated  that  the  end 
faces  tilted  instead  of  remaining  perpendicular  to  the  undeformed 
centroidal  axis  as  the  B.E.  theory  requires.  The  total  burden  of 
meeting  the  requirement  of  aero  slope  at  the  displaced  end  was 
put  on  the  QUAD4  elements  which  formed  the  side  panels  of  the 
steel  tube.  That  is;  the  open  ended  tube  had  two  surfaces  that 
could  carry  such  bending  and  two  surfaces  unable  to  carry 
in-plane  shear  about  their  normals.  Even  those  that  picked  up 
such  bending  couldn't  transmit  this  moment  to  the  QUAD'S  on  the 
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perpendicular  surface,  so  an  inadequate  moment  developed  to 
produce  a  net  slope  of  zero  at  the  centroid.  By  letting  the  end 
points  displace  axially  they  were  in  no  position  to  create 
couples  to  satisfy  the  moments  for  zero  slope.  That  is  why  the 
model  with  the  B.E.  constraints  produced  inadequately  stiff 
sectional  properties  in  bending  and  torsion. 

Going  back  to  the  model  with  fully  constrained  ends,  the 
explanation  as  to  why  this  model  was  also  inadequate  for  simulat¬ 
ing  an  equivalent  beam  was  this.  Even  though  it  did  develop 
couples  which  formed  the  resisting  moment  for  zero  end  slope  by 
holding  the  axial  displacements  to  zero;  it  still  felt  the  defi¬ 
ciency  of  moments  about  the  normals  of  side  panels.  In  effect 
membrane  action  on  corner  displacements  alone  was  not  sufficient 
to  represent  the  true  structure  without  the  help  of  the  existing 
—  but  unrepresented  --  in  plane  shear  from  moments  about  the 
normals  of  the  panels. 

In  the  case  of  torsion  the  fully  ;onstrained  model  was 
invalid  because  it  developed  local  equilibrium  at  the  end  under¬ 
going  unit  rotation.  The  unit  rotation  about  the  axial  direction 
for  every  end  grid  point  was  inhibited  by  the  translational 
d.o.f.'s  being  held  to  zero.  The  deformation  became  a  scalloped 
pattern  instead  of  a  uniformly  rotated  face.  Representation  of 
torsion  with  the  B.E.  model  was  also  inadequate  because  it  re¬ 
quired,  but  didn't  get,  the  assistance  of  the  panels  on  all  four 
sides  to  carry  the  rotation  about  their  normals. 

Does  this  mean  that  if  no  attempt  were  made  to  model  the 
mast  as  an  equivalent  beam,  but  a  full  3-D  model  were  used,  that 
the  3-D  model  would  be  invalid?  Not  at  all.  What  it  shows  is 
that  the  3-D  model  is  ineffective  in  trying  to  conform  to  the 
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requirements  of  an  equivalent  beam  representation.  If  a  full  3-D 
plate  model  were  used  in  the  complete  representation  of  the  crane 
structure,  good  results  would  be  obtained. 

Since  the  attempt  is  to  economize  on  the  size  of  the 
model,  a  better  way  to  achieve  the  same  results  is  to  use  sub¬ 
structuring  and  condense  the  mast  to  equivalent  end  boundary  and 
intermediate  mass  points. 

The  spirit  in  which  this  paper  is  presented  is  to  publish 
failures  as  well  as  successes  to  help  analysts  avoid  retracing 
the  ground  that  has  already  been  plowed. 
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ABSTRACT 

A  nonlinear  elastic  force-displacement  relationship  is  used  to  calculate  the 
transient  impact  force  and  local  deformation  at  the  point  of  contact  between  impactor  and 
target.  The  nonlinear  analysis  and  transfer  function  capabilities  of  NASTRAN  are  used  to 
define  a  finite  element  model  that  behaves  globally  linearly  elastic,  and  locally  nonlinear 
elastic  to  model  the  local  contact  behavior. 

Results  are  presented  for  two  different  structures:  a  uniform  cylindrical  rod 
impacted  longitudinally;  and  an  orthotropic  plate  impacted  transversely.  Calculated 
impact  force  and  transient  structural  response  of  the  targets  are  shown  to  compare  well 
with  results  measured  in  experimental  tests. 


INTRODUCTION 

Aerospace  structures  are  subjected  to  impact  loading  from  a  variety  of  sources, 
including  dropped  tools,  runway  debris,  and  munitions.  In  some  advanced  materials,  even 
low  velocity  impact  can  cause  significant  structural  damage.  Therefore,  development  of 
accurate  means  of  calculating  structural  response  due  to  impact  loading  can  be  of  critical 
importance.  In  this  paper,  a  computational  technique  is  developed  to  predict  the  dynamic 
response  of  a  structure  to  low  velocity  elastic  impact. 

Structural  damage  due  to  impact  invariably  initiates  in  the  immediate  vicinity 
of  the  impact.  Therefore,  it  is  important  that  the  local  stress  field  in  the  region  of  contact 
be  calculated  accurately.  Hertz  [1]  derived  an  elasticity-based  force-displacement 
relationship  that  describes  contact  between  two  elastic  bodies.  The  Hertzian  contact  law  is 
given  by: 


F  =  Ko" 


(1) 


where 


F  =  elastic  contact  force 
K  =  contact  stiffness 
n  =  exponent 
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a  =  relative  displacement  (indentation)  between  impactor  and  target 
=  Ui  -  ut  (i  =  impactor,  t  =  target) 

The  exponent  n  was  shown  in  reference  [1]  to  have  the  value  of  3/2-  In  dynamic 
applications  such  as  this;  F,  u,  and  a  are  all  time— varying. 

During  low  velocity  impact,  where  impact  damage  is  confined  to  the  area 
immediately  around  the  point  of  contact,  areas  of  the  structure  remote  from  the  impact 
may  still  deform  in  a  linear  elastic  manner.  An  efficient  finite  element  model,  therefore, 
would  combine  a  linear  elastic  model  of  the  global  structure  with  a  non— linearly  elastic 
behavior  at  the  point  of  contact.  The  nonlinear  force-displacement  relationship  in 
equation  (l)  is  incorporated  into  a  linear  elastic  finite  element  model  (MSC/NASTRAN 
transient  solution  27,  COSMIC/NASTRAN  transient  solution  9)  by  using  a  NASTRAN 
transfer  function  definition  and  nonlinear  analysis  capability.  In  the  following  section,  the 
Hertz  contact  law  is  discussed  in  addition  to  a  method  of  incorporating  it  into  NASTRAN. 
Impact  loading  of  two  different  structures  is  then  analyzed.  The  first  problem  is  a 
one-dimensional  rod  of  uniform  cross  section  impacted  longitudinally.  The  second  is  an 
orthotropic  plate  under  transverse  impact. 


CONTACT  LAW 

In  reference  [2]  Hertz  derived  the  force-displacement  relationship  for  two 
spherical  isotropic  elastic  bodies  of  radius  rj,  and  r2  in  contact: 

F  =  K  a3/2  (2) 


where 


is  the  contact  stiffness  and 


4 

ri  r2  ki  k2 

~  ~r 

r  i  +  r2  k  j  +  k2 

* 

E. 

k-  =  - 

C-) 

T— < 

II 

•t* 

r* 

1  -U1. 

(3) 


(4) 


where  E-  and  u-  are  the  Young’s  modulus  and  poisson’s  ratio,  respectively,  and  the 

subscripts  1  and  2  refer  to  each  of  the  spheres.  When  a  spherical  impactor  contacts  a  flat 
target,  (3)  simplifies  to 


kj  kt 
ki  +  kt 


(5) 


where  i  and  t  represent  the  impactor  and  target  respectively  and  the  kt  and  k\  are  given  by: 
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kt - S-J-  (6) 

1  -I't 

ki  =  Ei  i  (7) 

1  -I/i 

In  equation  (2),  a  is  the  local  indentation  at  the  contact  point,  shown 
schematically  in  figure  1.  We  have: 


a  =  ui-ut  (8) 

where  a  is  the  relative  local  displacement  between  impactor  and  target  at  the  point  of 
contact. 


NASTRAN  Implementation 

The  non-linear  local  behavior  was  incorporated  into  the  NASTRAN  finite 
element  model  as  follows: 

The  impactor  is  modeled  as  a  lumped  mass  just  touching  the  target  at  t=0  and 
with  an  initial  velocity  towards  the  target.  The  difference  between  the  displacement  of  this 
lumped  mass  and  the  displacement  of  the  target  is  the  indentation,  a.  The  modeling  of  the 
contact  between  impactor  and  target  is  performed  by  utilizing  the  transfer  function  card, 
TF,  and  the  nonlinear  force  card,  NOLIN3.  The  TF  card  acts  as  a  dynamic  multipoint 
constraint,  relating  the  displacement,  velocity  and  acceleration  of  several  independent 
degrees  of  freedom  to  a  dependant  degree  of  freedom.  In  the  work  discussed  here,  only 
displacement  relationships  were  used.  On  the  TF  card  coefficients  of  the  following 
equation  are  specified  [3]. 


(B0  +  B,p  +  B2p2)udep  +  X(AJ0+  AJlP  +  AJ2p2)uJind  =  0  (9) 

j  =  i 

where 


B0,  Bu  B2  =  the  coefficients  for  the  dependant  degree  of  freedom 

Aq,  Ap  A2  =  the  coefficients  for  the  independent  degrees  of  freedom 
udep  =  displacement  of  the  dependant  degree  of  freedom 

uind  =  displacements  of  the  independent  degrees  of  freedom 
n  =  the  number  of  independent  degrees  of  freedom 
p  =  the  differential  operator  — ,  and  p2  = 

For  this  analysis,  the  equation  would  appear: 


S 
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(l.O)llextra  point  +  J^(—1-O)uimpactor  +  (l.O)utarget 


=  0 


(10) 


that  is 


n 

®1>  ®2»  ^2 


B 

A 


=  2 


An  = 


0.0  (j  =  1,  n) 

1.0 

-1.0 

1.0 


The  resulting  equation  defines  the  indentation  at  every  time  step  and  assigns  the  value  to 
an  EPOINT.  The  EPOINT,  or  extra  point,  is  used  as  a  nonstructural  variable  that 
provides  a  place  to  store  the  value  of  the  indentation.  The  EPOINT  is  provided  as  input 
to  the  NOLIN3  card. 


The  NOLIN3  card  is  the  means  of  applying  the  time-dependent  nonlinear  load 
based  on  the  indentation.  The  NOLIN3  card  has  the  form: 


where 


[  S(x(t))A,  x(t)  >  0 


P(t) 


P(t)  =  is 
S  =  is 
x(t)  =  is 
A  =  is 


0  ,  x(t)  <  0 


the  resulting  nonlinear  force 
a  scale  factor 

the  displacement  or  velocity  of  a  degree  of  freedom 
an  amplification  factor 


(11) 


In  modeling  of  the  impact,  we  define  x(t)  to  be  the  displacement  of  the  EPOINT,  S  to  be 
the  Hertzian  stiffness,  and  A  to  be  3/ 2,  as  given  in  equation  (2).  Recall  that  the 
displacement  of  the  EPOINT  is  really  the  indentation  as  obtained  from  the  TF  card.  The 
resulting  function  then  has  the  form: 


p(t>  = 


K(0(t))3/’,  o(t)  >  0 
0  ,  a(t)  <  0 


(12) 
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Note  that  when  a  is  less  than  or  equal  to  zero  (ie.  the  target  and  the  impactor  are  out  of 
contact)  then  the  force  is  also  zero.  Two  NOLIN3  cards  are  used,  one  to  apply  the  impact 
force  to  the  target  and  the  other  to  apply  the  same  force  to  the  impactor  in  the  opposite 
direction  of  its  initial  velocity.  This  methodology  allows  the  impactor  to  slow  with 
increasing  impact  force  and  eventually  to  unload  the  target  as  the  impactor  begins  to  travel 
in  the  opposite  direction,  away  from  the  target. 


RESULTS 

One  Dimensional  Rod 

The  first  problem  analyzed  is  the  longitudinal  impact  of  a  steel  ball  on  a  long 
aluminum  rod  of  constant  cross  section.  Geometry  and  material  properties  of  the  impactor 
and  target  are  given  in  figure  1.  The  problem  was  modeled  using  144  1-D  rod  elements 
with  each  grid  point  having  a  single  longitudinal  degree  of  freedom.  Two  more  degrees  of 
freedom  were  used  to  model  the  impact  dynamics,  resulting  in  a  total  of  147  degrees  of 
freedom.  A  single  lumped  mass  with  an  initial  velocity  was  used  to  represent  the  impactor. 
The  Hertzian  force-displacement  relationship  in  equation  (1)  was  prescribed  using  the 
NASTRAN  NOLIN3  card,  as  shown  in  the  example  input  deck  in  the  appendix. 

The  calculated  impact  force  history  compares  well  with  experimentally 
determined  values  [4],  as  shown  in  figure  2.  The  calculated  strain  response  at  the  midpoint 
of  the  target  bar  is  compared  with  measured  values  in  figure  3.  The  sign  reversal  of  the 
second  pulse  is  caused  by  the  reflected  tensile  stress  wave  generated  by  the  incident 
compressive  wave  reaching  the  free  end  of  the  bar  [5]. 

Some  insight  into  the  timing  and  the  location  of  the  impact-induced  structural 
failure  can  be  gained  by  tracking  the  distribution  of  energy  in  the  impactor  and  the  target, 


as  shown  in  figure  4.  The  energy  balance  can  be  expressed  as: 

Ut  =  KEi  +  SEi  +  KEt  4  3Et  (13) 

where 

Ut  =  total  energy  in  system 

2 

KEi  =  impactor  kinetic  energy  =  l/2  mvj  (14) 

SEi  =  impactor  strain  energy  =  J F(a)  da  =  2/s  Ka  (15) 

KEt  =  target  kinetic  energy  (16) 

n  - 1 

SEt  =  strain  energy  of  target  (17) 


1  2 

-j—  kj  (  6j  —  )  (n  =  number  of  elements) 


The  total  energy  in  the  system,  Ut,  is  divided  between  the  kinetic  energy  and 
the  strain  energy  of  the  target  and  impactor  in  a  time-varying  manner.  Because  damping 
effects  are  not  considered,  the  total  system  energy  is  constant  and  equal  to  the  initial 
kinetic  energy  of  the  impactor.  The  strain  energy  of  the  impactor  is  non-zero  only  during 
the  contact  interval  (0  <  tC/L  <  0.4,  where  t  =  time,  L  =  the  length  of  the  bar,  and  C  = 
the  wave  speed  in  the  bar)  and  peaks  when  the  contact  force  is  greatest,  approximately 
halfway  through  the  contact  interval.  The  kinetic  energy  of  the  impactor  decreases  rapidly 
as  the  impactor  slows  during  contact  with  the  target.  Eventually,  at  tC/L  =  0.25,  the 
impactor  velocity  (and  therefore  its  kinetic  energy)  decreases  to  zero  and  the  elastic 
rebound  begins.  The  kinetic  energy  of  the  impactor  never  returns  to  its  initial  level 
because  approximately  80%  of  the  energy  has  been  transferred  to  the  target  in  the  form  of 
strain  energy  and  kinetic  energy.  The  strain  and  kinetic  energies  in  the  target  both 
increase  rapidly  during  the  contact  with  the  impactor  and  remain  constant  after  contact 
has  ended  (tC/L  >  0.4).  Both  strain  and  kinetic  energies  maintain  equal  and  constant 
values  until  the  compressive  stress  wave  generated  by  the  impact  reaches  the  far  end  of  the 
free-free  bar  (tC/L  =  1.0).  A  tensile  stress  wave  is  generated  when  the  compressive  pulse 
reflects  from  the  stress  free  boundary  [5].  The  superposition  of  the  incident  and  reflected 
pulses  momentarily  leaves  the  bar  stress-free  which  causes  the  strain  energy  to  decrease  to 
zero.  The  kinetic  energy  simultaneously  increases,  maintaining  a  conservation  of  total 
energy.  The  reflection  process  is  repeated  at  tC/L  =  2.0,  when  the  reflected  pulse  returns 
to  the  other  end  of  the  bar.  Similar  energy  dissipation  diagrams  may  prove  useful  in 
analyzing  dynamic  failure  of  more  complex  structures. 


Composite  Plate 

The  low  velocity  transverse  impact  of  a  composite  plate  made  from  Scotchply 
1003  prepreg  [6]  is  now  analyzed.  The  problem  is  depicted  schematically  in  figure  5,  and  is 
described  in  detail  in  references  [7,8].  A  modified  Hertzian  contact  stiffness  has  been 
proposed  [9]  for  application  to  composite  materials.  Specifically,  equation  (6)  is  replaced 
by 


kt  =  E33t  (18) 

where  E33t  is  the  transverse  modulus  of  the  plate.  Plate  membrane  and  bending  stiffness 
material  properties  were  calculated  using  the  COBSTRAN  (Composite  Blade  Structural 
Analyzer)  computer  code  [10]  which  calculates  elastic  moduli  of  composite  materials  from 
known  constituent  properties  and  laminate  ply  orientations. 

A  uniform  square  mesh  of  QUAD4  elements  was  used  to  model  the  15.24  cm  * 
15.24  cm  (6  in  *  6  in)  target  plate.  A  mesh  convergence  study  was  performed  to  establish 
the  degree  of  mesh  refinement  necessary  to  arrive  at  a  numerically  converged  solution. 
Three  different  meshes  were  considered,  25  *  25,  49  x  49,  and  61  *  61  elements.  Of  these, 
the  latter  two  produced  essentially  the  same  strain  response  for  a  given  impact  velocity  and 
were  therefore  considered  to  be  converged  solutions.  The  results  presented  here  were 
therefore  calculated  using  the  49  *  49  element  model.  Five  degrees  of  freedom  (ux,  uy,  uz, 
0X  and  0y)  were  used  at  each  nodal  point,  giving  the  model  a  total  of  11510  degrees  of 
freedom.  The  problem  was  solved  on  a  Cray  XMP  in  52  CPU  minutes. 
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The  impactor  used  in  the  tests  [7.8]  was  a  uniform  2.54  cm  (1  in)  long, 
blunt-ended  steel  rod  of  radius  0.047625  cm  (3/i6  in).  In  the  analysis  a  contact  radius  of 
0.047625  cm  (3/i6  in)  was  assumed  in  the  Hertzian  contact  stiffness  calculations.  The 
calculated  impact  force  history  is  shown  in  figure  6.  Although  no  direct  measurement  of 
the  impact  force  was  obtained  experimentally,  the  contact  time  was  measured  [8]  and  found 
to  be  204  microseconds.  This  is  in  good  agreement  with  the  calculated  result.  Figure  6 
also  shows  that  a  secondary  impact  occurs  during  the  latter  half  of  the  contact  interval  (t 
=  175  fise c),  probably  due  to  the  vibration  of  the  target  plate  during  contact  with  the 
projectile. 


The  resulting  displacement  response  of  the  plate  is  shown  in  figure  7,  where  it 
has  been  assumed  that  no  damage  occurs  in  the  target  during  contact  with  the  impactor. 
This  assumption  is  valid  based  on  the  available  test  data.  Ultrasonic  C-scans  of  the 
specimens  after  impact  indicate  that  this  level  of  impactor  kinetic  energy  (10  Joules)  is 
very  near  the  threshold  energy  level  required  to  cause  damage  [8]  in  specimens  of  this 
layup.  As  a  result,  very  little  damage  occurs  at  this  impactor  velocity.  The  anisotropic 
bending  stiffness  of  the  target  (figure  5)  is  evident  from  the  elliptical  displacement 
contours,  as  the  flexural  disturbance  travels  faster  in  the  stiffer  direction  (figure  7). 

The  strain  response  at  gage  A  is  compared  to  the  calculated  response  in  figure  8. 
The  two  curves  are  similar  in  amplitude  and  duration  but  the  calculated  strain  appears  to 
lag  the  measured  values  by  approximately  25  microseconds.  This  may  be  due  to  the 
difficulty  in  establishing  experimentally  the  precise  time  at  which  contact  occurs  based  on 
strain  gage  readings  taken  at  some  distance  from  the  point  of  contact.  The  comparison 
shown  in  figure  9  for  gage  B  likewise  shows  a  time  shift  of  approximately  25  microseconds 
between  the  measured  and  the  calculated  response.  The  amplitude  and  duration  of  the 
calculated  strain  response  correlate  quite  well  with  the  measured  signal. 


SUMMARY 

A  simple  means  of  modeling  low  velocity,  non-damaging  impact  using 
NASTRAN  was  demonstrated.  A  nonlinear  elastic  contact  model  was  included  in  the  finite 
element  analysis  using  NASTRAN  transfer  function  definitions  and  nonlinear  analysis 
capabilities.  The  same  contact  law  was  used  to  define  the  force-indentation  relationship 
for  two  different  impactor/target  combinations.  Results  in  both  cases  showed  that  the 
impact  force  and  resulting  transient  structural  response  of  the  target  compared  well  with 
experimentally  measured  values. 
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Figure  2:  Impact  Force  for  Bar 
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Figure  3: 


Specimen 
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Figure  5:  Composite  Plate  Impact 
Specimen  Configuration 
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Figure  7:  Calculated  Displacement  Response  for  Composite  Plate 
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Figure  8:  Strain  Response 
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DAN  TROWBRIDGE 


ANALEX  -  STRUCTRAL  MECHANICS  BRANCH 


ID  TRANS, LOAD 
APP  DISP 
TIME  60 
SOL  9 
CENO 

TITLE  ■  COSMIC:  TRANSIENT  RESPONSE  ANALYSIS:  HERTZIAN  IMPACT  FF 
SUBTITLE  -  36"  AL.  ROD  5/8  STEEL  BALL  V0-S2.1  IN/S 

LABEL  «  ROD  . - . — ■ - «—=—=■=*■=  •  < -  IMPACT 

i  NONLINEAR  LOAD 

NONLINEAR  «  5 

S  INITIAL  CONDITIONS  SET 

IC  -  1 

TFL-111 

SPC  -  4 

TSTEP  -  7 

S  OUTPUT  STUFF 

SET  30  -  1.72.73.999,1001 
NLLOAD  -  30 
STRESS (PR I NT)  -  30 
DISP(PRINT)  -  30 
BEGIN  BULK 

$ . . . * . . . . . . . 

S  EXTRA  POINT  «=  INDENTATION 

EPOINT, 1001 

GRID. 999. .-0.3125,0.0,0.0 
GRID.1 . .0.0, 0.0, 0.0 
-0*4).*(1).-.*(«.25).— 

CROD,  1.  1.  1.  2 
-(H3)..(1),-,.(1)..(1) 

$  LUMP  MASS  OF  IMP ACTOR 

CONM2. 200. 999. 0.9. 587-5, 0.0, 0.0, 0.0. ,  +CON2-2 

+CON2-2 , 3 . 745-6 , .3.745-6, . .3.745-6 

S  MATERIAL  PROPERTIES 

PROO.1 .11.0.196.6.14-3.0.25 

MAT1 ,11 ,10.0+6, .0.33.2.5-4, . . .+MAT1-1 

+MAT1-1 , 35 . 0E6 , 36 . 0E6 , 27 . 0E6 

$  BOUNDRY  CONDITIONS 

SPC 1 . 4 . 23456 . 1 . THRU .145 

$  REMOVE  DEGREES  OF  FREEDOM  FROM  IMPACTOR 

5 PCI  4  23456  999 

$  '  TRANSFER  FUNCTION  TO  DEFINE  INDENTATION 

TF.111 .1001 .0.+1 .0.0. 0.0.0, . ,+TF-l 

+TF-1 .999.1 .-1 .0.0. 0.0.0, , . .+TF-2 

+TF— 2 .1,1, 1.0, 0.0, 0.0 

*  TIMING 

TSTEP. 7. 2500. 2. 0-7, 25 

$  LOAD  DEPENDENT  ON  DISPLACEMENT  OF  IMPACTOR 

NOLIN3,  5.  1.  1.  6.24+6,  1001.  1.  1.5 
$  SLOW  DOWN  IMPACTOR 

NOLIN3.  5.  999.  1.  -6.24+6.  1001,  1.  1.5 

$  INITIAL  CONDITIONS:  IMPACTOR  VELOCITY  -  62.1  IN/SEC 

TIC. 1 .999,1 .0.0.62.1 

ENDDATA 
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ANAIEX  -  STRUCTRAl  MECHANICS  BRANCH 


ID  IMPACT. PLATE 
APP  OISP 
TIME  120 
SOL  27 
CEND 

TITLE  -  IMPACT  OF  PLATE  49X49  :  CENTERED  ELEMENT 
SU8TITLE  -  TRANSIENT  ANALYSIS:  FIXED-FIXED:  NO  SYMMETRY 

MrU  *  I 

IC  -  3 

NONLINEAR  -  5 
TSTEP  -  1 
TFL  -  111 

SET  15  -  999.2525.2526.2625.2626 
NLLOAD  -  15 

SET  20  -  2525.3325.4125 
STRESS  -  20 
BEGIN  BULK 

$  •*••  EXTRA  POINT  TO  HOLD  INDENTATION  •••••*... . .  .  .  . 

EPOINT, 10001  ****•• 

$  •*••  IMPACTOR  •*.  3/8  IN  DIAMETER  . . . 

GRIO. 999, .0.0,0. 0,-0. 1875 

CONM2 , 200, 999 , 0 , 8 . 096— 5 , 0 . 0 , 0 . 0  0  0  +CON2— 2 

+CON2-2. 7. 459-6,. 7. 459-6... 1.423-6  ' 

o  *  *  *  *  *  *  ?R1DS  AN0  CQUAtH  ELEMENTS  OEFINING  THE  PLATE  GO  HERE 

PSHELL.I.ID^^Te^1^-*-  mT2  CARD$  °ENERATED  BY  C0BSTRAN 

^iTii]I!8E-iitl69i4it?50Eli5E~03•2*8E+e6■_3-<E■02'5-7E+05'1-8E-0+•+A10, 

MAT2.201 . 5 . 7E+06 ,2.9 E+05 , — 1 .9E-04.1 .4E+06 .-3.8E-03.5  7E+05 
5  BOUNDRY  CONDITIONS 

SPC1.  1.  123456.  101,  THRU,  150 
SPC1.  1,  123456.  5001,  THRU,  505e 
SPC1,  1,  123456.  101 
-.=.-.•100 
**48 

SPC1 ,  1.  123456.  150 
=.-.-,• 100 
-48 

SPC1.  1.  12456.  999 
GRDSET . 6 

5  TIME  STEP  INFO 

TSTEP.  1.2000.  1 .0—7,  10 

WLIW.  5.2525.  50“.«M?“Se".T"5E  °ISPUCO,£W  °F  ‘“’“TOR 

NOLIN3.  5.2526.  3. +1.945+5.  10001,  0,  1.5 

NOLIN3,  5.2625.  3. +1.945+5.  10001,  0,  1.5 

NOLIN3.  5.2626,  3. +1.945+5.  10001,  0,  1.5 

$  SLOW  DOWN  IMPACTOR 

N0LIN3.  5.  999,  3.-7.779+5,  10001.  0,  1.5 

TC  111  «  TRANSFER  FUNCTION  TO  CALCULATE  INDENTATION 

TF. 1 1 1 , 10001 ,0.+l .0.00.0,00.0, , ,+TF-l 

+TF-1 .999.3.-1.0.00.0.00.0, , , , +TF-2 

+TF-2. 2525, 3. +0.25. 00. 0,00.0, , , .+TF-3 

+TF-3. 2526. 3. +0.25. 00. 0,00.0, , , .+TF-4 

+TF-4, 2625. 3. +0.25. 00. 0,00.0, , , .+TF-5 

+TF-5. 2626, 3. +0.25, 00. 0,00.0 

nC.3.999.3.0!!!!u70L  <;W0UlWS:  IMPACT0R  V*-OC!TY  -  1470  IN/SEC  (122.5  FT/SEC) 
ENODATA 


134 


TRANSITIONING  OF  POWER  FLOW  IN 
BEAM  MODELS  WITH  BENDS 


i>y 

Stephen  A.  Hambric 
Applied  Mathematics  Division  (184) 
David  Taylor  Research  Center 
Betliesdn,  MD  20084-5000 


ABSTRACT 

The  propagation  of  power  flow  through  a  dynamically 
loaded  beam  model  with  90  degree  bends  is  investigated  using 
NASTRAN  and  McPOW.  The  transitioning  of  power  flow 
types  (axial,  torsional,  and  flexural)  is  observed  throughout 
the  structure.  To  get  accurate  calculations  of  the  torsional 
response  of  beams  using  NASTRAN,  torsional  inertia  effects 
had  to  be  added  to  the  mass  matrix  calculation  section  of  the 
program.  Also,  mass  effects  were  included  in  the  calculation 
of  BAR  forces  to  improve  the  continuity  of  power  flow 
between  elements.  The  importance  of  including  all  types  of 
power  flow  in  an  analysis,  rather  than  only  flexural  power,  is 
indicated  by  the  example.  Trying  to  interpret  power  flow 
results  that  only  consider  flexural  components  in  even  a 
moderately  complex  problem  will  result  in  incorrect 
conclusions  concerning  the  total  power  flow  field. 

INTRODUCTION 

Methods  for  calculating  power  flows  in  dynamically  loaded  finite 
element  models  using  NASTRAN  (Rigid  Format  8  -  Direct  Frequency 
Response)  and  McPOW  (Mechanical  POWer)  were  developed  previously.1 
The  power  flow  equations  for  beam  elements  derived  in  that  paper  included  all 
forms  of  dynamic  energy  propagation:  flexural,  longitudinal  (or  axial),  and 
torsional.  The  flexural  waves  were  split  into  shear  and  moment  components. 

The  majority  of  procedures  employed  in  other  studies  (sec  the  list  of 
references  in  Ilambric1)  only  consider  flexural  vibration  in  their  calculations  of 
power  flow.  This  can  be  dangerous  if  an  analyst  is  investigating  the  energy 
propagation  characteristics  of  a  complex  structure.  Though  flexural  vibration 
is  in  most  cases  the  dominant  response  in  a  dynamically  excited  beam, 
different  kinds  of  propagation  will  occur  in  structures  with  even  a  small  degree 
of  complexity,  such  as  a  simple  beam  model  with  90-degree  bends. 
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Such  a  model  is  tested  here  using  a  frequency  range  spanning  several 
resonances  and  types  of  motion.  Plots  showing  the  contributions  of  the 
different  forms  of  power  llow  to  the  total  power  travelling  through  the  system 
are  shown,  and  illustrate  the  importance  of  all  types  of  energy  propagation  to 
the  power  flow  method. 

To  improve  the  accuracy  of  both  the  finite  element  solution  and  the 
power  flow  solution  of  the  problem,  a  few  modifications  were  made  to 
NASTRAN  and  McPOW.  First,  to  show  the  importance  of  torsional  power 
flow,  a  capability  to  calculate  dynamic  torsional  forces  and  corresponding 
angular  velocities  is  required.  Therefore,  torsional  inertias  were  added  to  the 
coupled  mass  matrix  formulation  of  the  BAR  element.  Also,  since  the  beam 
element  force  calculation  algorithm  in  NASTRAN  considers  only  stiffness 
effects,  mass  and  damping  effects  had  to  be  added  to  McPOW  to  modify  the 
element  forces. 

METHODOLOGY 

The  procedure  for  solving  for  the  power  llow  field  in  a  finite  element 
model  using  NASTRAN  and  McPOW  is: 

1.  Run  Rigid  Format  8  (Direct  Frequency  Response)  on  a  NASTRAN  data 
deck  (using  the  ALTER  statements  shown  in  Ref.  1  to  output  force  and 
velocity  data  blocks  to  the  OUTPUT2  file).  Coupled  mass  formulations 
should  always  be  used. 

2.  Run  McPOW  using  the  binary  data  in  the  OUTPUT2  file  as  input. 

General  Methods 

A  typical  power  flow  cycle  is  shown  in  Fig.  1.  The  figure  shows  an 
arbitrary  structure  mounted  to  a  connecting  structure  by  a  spring  and  damper 
coupling.  A  dynamic  load  is  applied,  and  energy  Jlows  into  the  structure  at  the 
load  point.  The  input  power  then  Hows  through  the  structure  along  multiple 
llow  paths  denoted  by  arrows  whose  lengths  represent  power  llow  magnitudes. 
As  the  energy  Hows  toward  the  mounting,  it  is  dissipated  by  material  damping 
and  sound  radiation  into  a  surrounding  medium,  and  the  llow  arrows  shorten. 
The  flow  and  dissipation  processes  continue  until  the  remaining  energy  exits 
the  structure  through  the  mounting  and  Hows  into  the  connecting  structure. 
Though  only  one  power  entry  point  and  one  exit  point  are  shown  in  this 
drawing,  multiple  loads  and  mountings  may  exist.  A  classic  text  which 
describes  the  flow  of  structure-borne  sound  is  the  book  by  Cremcr,  Hcckl,  and 
Ungar.2 

The  structural  dynamics  problem  may  be  solved  using  NASTRAN.  The 
structure  may  be  modeled  with  various  element  types;  mountings  are  modeled 
with  scalar  spring,  damping,  and  mass  elements.  Constraints  and  loads  arc 
directly  applied.  The  steady-slate  response  for  the  model  is  solved  for  a  given 
excitation  frequency,  and  the  power  llow  variables  arc  calculated. 
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Fig.  1.  Sample  Power  Flow  Diagram. 

Power  is  defined  as  the  time-averaged  product  of  a  force  with  the  in- 
phase  component  of  velocity  in  the  direction  of  the  force.  For  time-harmonic 
analysis,  where  complex  numbers  are  used,  this  calculation  may  be  visualized 
as  taking  the  dot  product  of  the  force  and  velocity  phasors.  (There  is  no  factor 
1/2  in  the  following  power  equations  if  the  assumption  that  forces  and 
velocities  are  "effective"  values  rather  than  amplitudes  is  made.  With  this 
assumption,  consistency  is  maintained,  and  there  is  no  mixing  of  effective  and 
peak  quantities  in  this  formulation.) 

Multiplying  one  complex  number  by  the  in-phase  part  of  another 
complex  number  is  the  same  operation  as  multiplying  the  first  number  by  the 
complex  conjugate  of  the  other  number  and  taking  the  real  part  of  the  result. 
Therefore  a  general  formula  for  power  How  in  a  structure  is 

Power  =  Real  [Fv*  ] ,  ( 1 ) 

where 

F  =  force,  and 

v  =  complex  conjugate  of  velocity. 
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Power  Flow  Equations 

The  equations  for  power  flows  in  BAR  elements  are  repeated  here.  A 
diagram  of  the  BAR  element  and  its  NASTRAN  force  output  conventions  is 
shown  in  Fig.  2,  where  Plane  1  is  vertical  and  Plane  2  is  horizontal. 


Fig.  2.  The  BAR  Element 

Since  a  beam  is  a  one-dimensional  element,  energy  flows  in  only  one  direction: 
in  the  local  x  direction,  or  along  the  length  of  the  beam.  The  total  power  flow 
for  a  beam  element  is 

Px  =  Real  [  —  (l’’x vx  + V i  vy  +  V 2 vz +T0X  —M2 by  +M 1 ) ] ,  (2) 

where 

Fx  =  axial  force, 

Vy  —  shear  force  in  y  direction, 

V2  =  shear  force  in  z  direction, 

T  =  torsion  about  x, 

M2  =  bending  moment  about  y, 

Mi  =  bending  moment  about  z, 

Vj  =  translational  velocities  in  direction  i,  and 
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0;  =  rotational  velocities  about  axis  i. 


The  negative  sign  on  the  result  comes  from  force  and  displacement  direction 
conventions  for  the  element.  The  negative  sign  on  the  Mo  term  reflects  the 
NASTRAN  force  output  convention.  In  Fig.  2,  Mo  is  shown  as  positive  in  the 
opposite  sense  to  0y.  Therefore,  M20y  is  opposite  in  sign  to  the  other  power 
flow  components. 

NASTRAN  Modifications 

Torsional  Inertias 

NASTRAN  currently  does  not  consider  torsional  inertias  in  its  beam 
element  formulation.  Therefore,  all  torsional  results  (angular  displacements 
and  torques)  are  based  on  stiffness  only,  and  are  essentially  those  of  a  static 
problem  solution.  To  remedy  this,  torsional  inertias  were  added  to  the 
coupled  mass  formulation.  At  the  point  in  NASTRAN  where  the  basic 
element  mass  matrix  is  formed,  no  consideration  is  given  to  beam  offsets  or 
beam  orientation;  all  mass  coefficients  (as  well  as  stiffness)  are  calculated  in 
the  local  beam  coordinate  system. 

The  torsional  mass  moment  of  inertia  of  a  beam  is  p  L  Jx/2,  where  p  is 
the  mass  density,  L  is  the  beam  length,  and  Jx  is  the  polar  area  moment  of 
inertia.  In  the  standard  consistent  mass  matrix  for  a  beam,13  this  value  is 
broken  up  into  2/3  and  1/3  components;  2/3  of  the  value  is  placed  at  the 
diagonal,  and  1/3  is  placed  at  the  coupled  degree  of  freedom  (the  node  on  the 
other  end  of  the  beam).  The  same  fractions  are  used  for  the  translational,  or 
axial  masses.  In  NASTRAN,  however,  the  coupled  mass  formulation  uses  an 
average  of  lumped  and  consistent  formulations  to  reduce  error.  This  average 
changes  the  components  to  5/6  and  1/6  of  the  total  value.  Since  these  values 
are  currently  used  for  the  axial  masses  in  NASTRAN,  they  were  also  used  for 
the  torsional  inertias. 

Element  Force  Calculations 

NASTRAN  element  forces  are  currently  calculated  by  multiplying 
element  stiffness  matrices  by  element  displacement  vectors.  Both  damping 
and  mass  effects  are  ignored.  The  damping  in  a  stiffness  element  is  actually  in 
the  form  of  a  loss  factor,  which  generates  a  complex  stiffness  matrix.  All 
stiffness  terms  are  multiplied  by  1.0  +  i //.  For  most  dynamic  analyses, 
neglecting  the  h]  term  is  acceptable  since  it  is  generally  small.  For  a  power 
flow  analysis  of  a  highly  reverberant  structure,  however,  ignoring  the  loss 
factor  is  disastrous.  In  a  highly  reverberant  structure,  the  force  and  velocity  at 
a  given  point  are  close  to  90  degrees  out  of  phase.  Since  power  flow  is  defined 
as  the  dot  product  of  these  two  components,  a  small  change  in  the  phase  of 
the  force  has  large  effects  on  the  calculated  element  powers. 

Neglecting  the  element  mass  matrices,  whose  components  arc  several 
orders  of  magnitude  less  than  those  of  the  stiffness  matrices,  has  less  drastic 
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effects  on  the  power  flow  solution,  since  at  low  frequencies  the  masses  will 
have  little  effect  on  the  force  calculations  (the  element  mass  matrix  is 
multiplied  by  —  oj  2  to  take  the  second  time  derivative  of  the  corresponding 
displacements).  However,  when  high  frequency  analyses  are  performed  on  a 
model,  the  — w  2  multiplying  factor  becomes  more  significant,  and  neglecting 
the  mass  contributions  will  cause  some  error  in  the  force  calculations.  Errors 
in  element  forces  cause  errors  in  element  power  flows. 

Including  these  missing  effects  in  NASTRAN  is  complicated  by  the  fact 
that  the  element  force  calculation  algorithm  splits  the  problem  into  real  and 
imaginary  parts.  The  element  stiffness  matrices  are  multiplied  by  the  real 
parts  of  the  displacement  vectors  to  calculate  real  force  components,  and  the 
process  is  repeated  for  the  imaginary  components.  Adding  an  imaginary  term 
to  the  stiffness  matrices  causes  new  terms  to  be  generated  in  the  multiplication 
(imaginary  stiffness  x  imaginary  displacement  and  imaginary  stiffness  x  real 
displacement).  There  is  also  no  frequency  dependence  in  the  current 
algorithm,  since  stiffness  are  frequency  independent.  Mass  matrices,  however, 
must  be  multiplied  by  the  — oj  2  term  mentioned  above,  so  they  must  be 
recalculated  for  every  frequency. 

To  avoid  these  complications,  the  element  force  calculations  were 
temporarily  moved  to  McPOW.  The  element  mass  and  complex  stiffness 
matrices  arc  recalculated  on  a  local  clement  level,  and  combined  with  local 
element  displacements  to  solve  for  element  forces.  A  force  vector  with  12 
entries  is  the  result;  shears  in  the  local  y  and  z  directions,  moments  about  the 
local  y  and  z  directions,  axial  forces,  and  torques  arc  solved  for  at  each  grid 
point.  In  NASTRAN  only  eight  forces  are  calculated,  because  only  moments 
are  calculated  at  both  ends  of  a  beam  element.  Beam  power  flows  are 
therefore  calculated  at  each  end  of  the  element  using  only  the  forces  at  that 
end  and  the  corresponding  grid  velocities.  The  average  of  the  powers  at  the 
ends  is  taken  to  find  an  element  power  llow. 

TEST  PROBLEM 


Problem  Statement 

The  beam  model  that  was  analyzed  is  shown  in  Fig.  3.  All  three 
sections  have  the  same  cross  section  and  material  properties.  Dashpots 
(DAMP2  elements)  of  value  106  were  applied  at  the  model’s  end  in  all  six 
degrees  of  freedom.  A  unit  load  was  applied  at  the  top  end  of  the  model  in 
the  longitudinal  direction  (along  the  -  z  axis)  over  a  frequency  range  of  1  to  250 
Hz  swept  in  1  Hz  increments.  The  finite  element  model  consists  of  152 
elements  and  153  grid  points.  Grid  and  clement  numbering  starts  at  the  left 
end  of  Link  1  and  proceeds  up  to  the  end  of  Link  3. 

Using  the  local  beam  clement  coordinate  systems  shown  in  Fig.  3,  the 
following  table  of  force  balances  at  the  corners  (Link  3  to  Link  2,  Link  2  to 
Link  1)  was  generated. 
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Fig.  3.  Tesi  Problem  Geometry 
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The  subscripts  on  the  shears  and  moments  refer  to  the  plane  in  which  the 
forces  occur  (sec  Fig.  2).  This  table  can  be  used  to  track  the  propagation  of 
power  flow  through  the  structure.  For  example,  the  longitudinal  power  input 
to  Link  3  will  travel  down  the  beam  in  axial  waves  to  the  first  bend  and 
become  shear  power  How  in  the  z  direction  in  Link  2.  This  shear  power  will 
interchange  with  moment  power  along  the  beam  (the  sum  of  the  shear  and 
moment  components  is  the  total  flexural  power  flow  in  the  beam).  Any  shear 
power  that  exists  at  the  lower  end  of  Link  2  will  transition  to  more  shear 
power  in  Link  1. 
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Results  and  Discussion 

The  computed  power  input  curve  over  the  excitation  frequency  range  is 
shown  in  Fig.  4.  The  power  input  peaks  correspond  to  various  resonances  in 
the  structure.  Most  are  flexural,  but  some  axial  and  torsional  modes  influence 
the  power  input  curve.  The  longitudinal  modes  of  Link  3  cause  power  input 
peaks  (at  190  Hz  and  above),  as  well  as  the  torsional  modes  of  Link  1  (at  151 
Hz  and  above). 

In  this  model,  the  power  flow  path  is  independent  of  frequency.  The 
total  power  must  always  flow  from  (lie  input  point  at  the  end  of  Link  3  to  the 
dampers  at  the  beginning  of  Link  1.  This  simplifies  the  interpretation  of  the 
results,  since  the  directions  of  total  power  flow  are  established. 

The  types  of  power  flow  in  a  given  link  are  not  so  well-defined. 
Whether  the  dominant  path  in  a  link  is  flexural,  axial,  or  torsional,  depends  on 
(he  motion  of  the  structure.  Fig.  5  shows  the  two  most  common  types  of 
motion  paths  for  this  problem.  The  displacement  field  of  Diagram  ]  occurs 
most  often.  The  axial  load  applied  to  Link  3  drives  the  entire  structure 
forward  and  backward  over  a  frequency  cycle.  The  dominant  power  flow  in 
Link  3  is  axial;  the  dominant  power  flow  in  Link  2  is  flexural;  and  torsional 
and  flexural  power  flows  are  dominant  in  Link  1,  since  the  input  load  applies 
both  a  torque  and  bending  moment  to  the  link. 

In  Diagram  2  of  Fig.  5  a  different  type  of  motion  is  shown.  The  axial 
load  still  drives  the  upper  half  of  the  structure  in  the  same  direction,  but  the 
lower  half  moves  in  the  opposite  direction.  This  type  of  motion  is  not  what 
one  would  expect  in  a  static  problem,  but  the  dynamic  characteristics  of  the 
structure  produce  this  type  of  motion  in  various  frequency  ranges. 

Due  to  this  motion  path,  the  axial  power  flow  travelling  down  Link  3 
becomes  flexural,  torsional,  and  axial  in  Link  2.  The  torsional  and  axial 
components  appear  because  the  link  is  twisted  and  stretched  by  the  opposite 
directions  of  motion  of  the  two  ends.  The  torsional  power  in  Link  2  becomes 
flexural  power  in  Plane  1  in  Link  1,  and  the  axial  power  in  Link  2  turns  into 
flexural  power  in  Plane  2  in  Link  1.  The  flexural  power  in  Link  2  becomes 
torsional  and  flexural  power  in  Link  1  as  before  (Diagram  1). 

Considering  these  modes  of  power  transitioning,  the  power  flow  plots  in 
Figs.  6-8  may  be  interpreted.  Each  plot  shows  the  contributions  of  flexural, 
torsional,  and  axial  power  flow  as  a  percentage  of  the  total  power  flow  in  the 
center  of  each  link. 

Fig.  8  shows  the  power  components  in  Link  3.  Since  the  input  power  is 
in  the  longitudinal  direction,  the  majority  of  the  power  in  this  link  is  axial.  At 
certain  frequencies,  the  percentage  of  axial  power  is  greater  than  100  percent. 
The  large  axial  percentage  arises  because  at  certain  frequencies,  reflected 
waves  carry  power  in  the  opposite  direction  (toward  the  load).  Three  flexural 
resonances  in  the  structure  cause  reflected  power  just  before  50  Hz,  along  with 
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Fig.  5.  Dominant  motion  paths  for  test  problem 


five  others  right  after  100  Hz.  Between  200  and  250  Hz,  some  flexural  and 
torsional  resonances  cause  more  reflected  powers. 

Fig.  7  shows  the  power  components  in  Link  2.  The  dominant  type  of 
power  is  the  flexural  component  in  Plane  1,  and  is  denoted  by  the  solid  curve. 
This  type  of  power  field  corresponds  to  the  motion  type  shown  in  Diagram  1  in 
Fig.  5.  However  at  certain  frequencies,  the  power  flow  pattern  of  Diagram  2 
becomes  dominant,  and  axial  and  torsional  power  become  important.  In  most 
cases,  the  axial  power  flows  forward  (away  from  the  load  point),  and  the 
torsional  power  is  backward  (reflected  toward  the  load  point).  These 
tendencies  occur  at  the  same  frequencies  as  the  reflected  power  waves  do  in 
Link  3  (shown  in  Fig.  8).  This  behavior  indicates  that  the  flexural  power  in 
Plane  1  and  the  torsional  power  cause  reflected  flexural  powers  in  Planes  2  and 
1  respectively  in  Link  3. 

Fig.  6  shows  the  power  distribution  in  Link  1.  In  this  case,  all  power 
components  are  positive,  implying  that  the  reflected  power  waves  in  Links  2 
and  3  originate  from  the  joint  connecting  Links  1  and  2.  In  Link  1,  flexure  in 
Plane  1  and  torsion  are  the  dominant  components  of  power  (low.  Flexural 
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Frequency 


motion  in  Plane  2  and  axial  motion  cause  power  peaks  at  the  same  frequencies 
observed  in  Figs.  7  and  8,  indicating  the  type  of  motion  shown  in  Diagram  2 
in  Fig.  5.  A  torsional  mode  in  Link  1  accounts  for  the  peak  in  the  torsion 
curve  at  150  Hz,  along  with  an  input  power  peak  at  the  same  frequency  (see 
Fig- 4). 

In  spite  of  the  large  variation  in  percentages  of  power  types  in  the  plots, 
all  the  power  curves  add  up  to  100  percent,  as  expected.  In  addition,  the  total 
power  ilow  in  the  structure  at  all  frequencies  is  at  a  maximum  at  the  load 
point,  and  smoothly  decreases  to  a  minimum  at  the  connection  point  to  the 
dampers.  The  steady  decrease  in  power  is  due  to  structural  damping.  The 
remaining  power  is  all  dissipated  by  the  connected  dampers. 

This  example  illustrates  the  importance  of  all  types  of  power 
components  in  a  power  flow  analysis.  Imagine  trying  to  discern  a  meaningful 
power  flow  field  from  only  flexural  powers  in  this  example.  The  detected 
powers  in  Link  3,  which  is  adjacent  to  the  input  load,  are  all  in  the  oposite 
direction,  or  toward  the  load.  In  Link  2,  the  analyst  would  see  a  sudden  jump 
in  power  to  values  that  are  higher  than  that  of  the  input  power.  Finally,  in 
Link  1,  sporadic  power  curves  with  values  near  the  input  power  at  frequencies 
below  100  Hz  and  values  near  zero  after  100  Hz  would  be  found.  Confusion 
would  surely  result,  with  erroneous  conclusions  soon  following.  Difficulties 
like  these  would  be  compounded  in  a  real  application  with  some  degree  of 
complexity. 


CONCLUSIONS 

The  modifications  made  to  NASTRAN  and  McPOW  are  critical  to  the 
powet  flow  method.  Without  torsional  inertias  applied  to  the  beam  element 
mass  matrices,  any  torsional  effects  in  a  dynamic  problem  are  sialic.  None  of 
the  torsional  power  flows  present  in  the  example  problem  would  exist,  causing 
incorrect  total  power  flow  fields.  Adding  mass  and  damping  effects  to  the 
element  force  calculation  algorithm  is  also  important.  In  a  reverberant 
structure  where  forces  and  velocities  are  nearly  90  degrees  out  of  phase  with 
each  other,  accurate  calculations  are  necessary  to  guarantee  good  power  flow 
results.  A  small  change  in  the  phase  of  an  element  force,  caused  by  neglecting 
the  material  loss  factor,  could  cause  large  errors  in  clement  power  flows. 
Also,  at  higher  frequencies,  element  mass  terms  can  become  significant  and 
affect  the  element  force  magnitudes,  and  hence  the  clement  power  magnitudes. 

The  addition  of  torsional  inertia  to  the  beam  element  mass  matrix 
formulation  was  straight-forward.  The  addition  of  damping  and  mass  effects  to 
the  element  force  calculation  routines,  however,  was  almost  impossible.  In 
fact,  the  changes  had  to  be  made  to  McPOW  instead  of  NASTRAN.  The 
implementation  difficulties  were  due  to  the  way  NASTRAN  handles  complex 
analysis:  the  solutions  are  broken  into  real  and  imaginary  parts.  When  the 
program  was  in  its  formative  stages,  UNIVAC  computers  were  supported. 
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The  UNIVAC,  unfortunately,  had  no  way  of  handling  double  precision 
complex  arithmetic.  Therefore,  no  complex  numbers  or  FORTRAN  complex 
functions  are  used  in  the  element  force  calculation  sections  of  the  program. 
With  this  approach,  a  simple  complex  calculation  like 
\—ur  [M]c  +  (1  4-  i7?)[K]e]  {d}c  must  be  split  up  into  four  calculations.  Also, 
since  (he  calculation  is  frequency-dependent,  the  NASTRAN  element  force 
subroutines  are  not  currently  able  to  handle  it.  Since  the  UNIVAC  has  all  but 
disappeared  from  the  COSMIC  NASTRAN  computing  arena  and  most 
modern  computers  support  double  precision  complex  arithmetic,  perhaps  the 
way  NASTRAN  handles  complex  problems  should  be  modified. 

The  importance  of  including  longitudinal  and  torsional  components  with 
flexural  ones  in  a  power  flow  analysis  was  shown  in  the  example  problem. 
Measuring  flexural  power  alone  will  not  give  an  accurate  indication  of  the  total 
power  flow  field  in  even  a  marginally  complex  problem.  In  the  case  of  the 
example  problem,  reflected  flexural  waves  actually  indicated  a  reversal  of 
power  flow  in  the  model,  where  the  direction  of  flexural  power  was  in  the 
opposite  direction  of  the  input  power.  Trying  to  interpret  power  flow  results 
that  only  consider  flexural  components  will  result  in  incorrect  conclusions 
concerning  the  total  power  flow  field. 
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ABSTRACT 

NASHUA  is  a  coupled  finite  element/boundary  element  capability  built  around 
NASTRAN  for  calculating  the  low  frequency  far-field  acoustic  pressure  field  radiated  or 
scattered  by  an  arbitrary,  submerged,  three-dimensional,  elastic  structure  subjected  to  either 
internal  time-harmonic  mechanical  loads  or  external  time-harmonic  incident  loadings.  This 
paper  describes  the  formulation  and  use  of  NASHUA  for  solving  such  structural  acoustics 
problems  when  the  structure  is  fluid-filled.  NASTRAN  is  used  to  generate  the  structural 
finite  element  model  and  to  perform  most  of  the  required  matrix  operations.  Both  fluid 
domains  are  modeled  using  the  boundary  element  capabililty  in  NASHUA,  whose  matrix 
formulation  (and  the  associated  NASTRAN  DMAP)  for  evacuated  structures  can  be  used 
with  suitable  interpretation  of  the  matrix  definitions.  After  computing  surface  pressures  and 
normal  velocities,  far-field  pressures  are  evaluated  using  an  asymptotic  form  of  the 
Helmholtz  exterior  integral  equation.  The  proposed  numerical  approach  is  validated  by 
comparing  the  acoustic  field  scattered  from  a  submerged  fluid-filled  spherical  thin  shell  to 
that  obtained  with  a  series  solution,  which  is  also  derived  in  this  paper. 


INTRODUCTION 

Two  basic  problems  in  computational  structural  acoustics  are  (1)  the  calculation  of  the 
acoustic  pressure  field  radiated  by  a  general  submerged  three-dimensional  elastic  structure 
subjected  to  internal  time-harmonic  loads,  and  (2)  the  calculation  of  the  acoustic  pressure 
scattered  by  such  a  structure  subjected  to  an  incident  time-harmonic  wavetrain.  The  most 
common,  as  well  as  the  most  accurate,  approach  for  solving  these  problems  at  low 
frequencies  is  to  couple  a  finite  element  model  of  the  structure  with  a  boundary  element 
model  of  the  surrounding  fluid.  This  is  the  approach  taken  by  NASHUA,  which  is  a 
boundary  clement  program  built  around  NASTRAN,  a  widely-used  finite  element  computer 
program  for  structural  dynamics. 

Several  previous  papers  (Ref.  1-4)  described  the  basic  formulation  and  development  for 
acoustic  radiation  and  scattering  from  evacuated  structures.  Here  we  describe  the 
formulation  and  use  of  NASHUA  for  modeling  submerged  structures  which  are  fluid-filled. 
Internal  fluid  can  occur  because  the  structure  is  free-flooded  or  contains  fluid-filled  tanks.  It 
is  possible  to  use  existing  NASTRAN  capability  to  model  the  interior  fluid  with  finite 
elements  (Ref.  5-7),  but  three-dimensional  models  with  large  numbers  of  fluid  degrees  of 
freedom  might  result.  An  attractive  alternative  to  the  fluid  finite  clement  model  is  to 
represent  the  contained  fluid  using  a  boundary  element  approach.  In  principle,  any  computer 
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program  capable  of  generating  the  appropriate  boundary  element  matrices  for  an  exterior 
fluid  is  also  capable  of  generating  such  matrices  for  the  complementary  region  (the  interior 
region).  NASTRAN’s  versatility  in  user-controlled  matrix  operations  (DMAP)  makes  the 
implementation  of  such  an  approach  straightforward. 


THEORETICAL  APPROACH 

The  basic  theoretical  development  for  NASHUA’S  radiation  and  scattering  approach 
for  evacuated  structures  has  been  presented  in  detail  previously  (Ref.  1-4).  For 
completeness,  this  paper  summarizes  that  approach  and  describes  the  changes  necessary  to 
model  the  interior  fluid  with  boundary  elements  in  the  same  procedure.  There  is  no 
requirement  that  the  interior  and  exterior  fluids  be  the  same. 


The  Surface  Solution  for  Evacuated  Structures 

Consider  any  submerged  three-dimensional,  evacuated  elastic  structure  subjected  to 
either  internal  time-harmonic  loads  or  an  external  time-harmonic  incident  wavetrain.  If  the 
structure  is  modeled  with  finite  elements  using  NASTRAN,  the  resulting  matrix  equation  of 
motion  can  be  written  as 

Zv  =  F  —  GAp,  (1) 

where  matrix  Z  (of  dimension  s  x  s)  is  the  structural  impedance,  vector  v  (s  x  r)  is  the 
complex  velocity  amplitude  for  all  structural  DOF  (wet  and  dry)  using  the  coordinate  systems 
selected  by  the  user,  vector  F  (s  x  r)  is  the  complex  amplitude  of  the  mechanical  forces 
applied  to  the  structure,  matrix  G  (s  x  f)  is  the  rectangular  transformation  of  direction 
cosines  to  transform  a  vector  of  outward  normal  forces  at  the  wet  points  to  a  vector  of 
forces  at  all  points  in  the  coordinate  systems  selected  by  the  user,  matrix  A  (f  x  f)  is  the 
diagonal  area  matrix  for  the  wet  surface,  and  vector  p  (f  x  r)  is  the  complex  amplitude  of 
total  pressures  (incident  +  scattered)  applied  at  the  wet  grid  points.  In  this  equation,  the 
time  dependence  exp(iwt)  has  been  suppressed.  In  the  above  dimensions,  s  denotes  the  total 
number  of  independent  structural  DOP  (wet  and  dry),  f  denotes  the  number  of  fluid  DOF 
(wet  points),  and  r  denotes  the  number  of  load  cases.  In  general,  the  surface  areas,  the 
normals,  and  the  transformation  matrix  G  are  obtained  in  NASHUA  from  the  NASTRAN 
calculation  of  the  load  vector  resulting  from  an  outwardly  directly  static  unit  pressure  load  on 
the  structure’s  wet  surface. 

In  Eq.  1,  the  structural  impedance  matrix  Z,  which  converts  velocity  to  force,  is  given 
by 

Z  =  (-arM  +  iwB  +  K)/(iw),  (2) 

where  M,  B,  and  K  are  the  structural  mass,  viscous  damping,  and  stiffness  matrices, 
respectively,  and  u  is  the  circular  frequency  of  excitation.  For  structures  with  a  nonzero  loss 
factor,  K  is  complex.  In  addition,  K  can  include  the  differential  stiffness  effects  of 
hydrostatic  pressure,  if  any  (Ref.  3).  A  standard  NASTRAN  finite  element  model  of  the 
structure  supplies  the  matrices  K,  M,  and  B. 

For  the  exterior  fluid  domain,  the  total  fluid  pressure  p  satisfies  the  Helmholtz 
differential  equation 
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V2p  +  k2p  =  0, 


(3) 


where  k  =  ulc  is  the  acoustic  wave  number,  and  c  is  the  fluid  sound  speed.  Equivalently,  for 
the  exterior  fluid,  p  is  the  solution  of  the  Helmholtz  integral  equations  (Ref.  8) 


fs  p(x)“^-dS  -  Js  q(x)D(r)dS  = 


p(x')/2-pi,  x'  on  S, 
p(x')-pi,  x'inE, 
-nr.  x'  in  I, 


(4) 


where  S,  E,  and  I  denote  the  surface,  exterior,  and  interior  domains,  respectively,  pi  is  the 
incident  free-field  pressure  (if  any),  r  is  the  distance  from  x  to  x'  (Fig.  1),  D  is  the  free-space 
Green’s  function 


D(r) 


q 


— iw/>vn. 


(5) 

(6) 


FLUID 


Pi 


Fig.  1.  Notation  for  Helmholtz  Integral  Equation 


p  is  the  fluid  mass  density,  and  vn  is  the  outward  normal  velocity  on  S.  As  shown  in  Fig.  1,  x 
in  Eq.  4  is  the  position  vector  for  a  typical  point  Pj  on  the  surface  S,  x'  is  the  position  vector 
for  the  point  P,  on  the  surface  or  in  the  exterior  field,  the  vector  r  =  x'  -  x,  and  n  is  the  unit 
outward  normal  at  Pj.  We  denote  the  lengths  of  the  vectors  x,  x',  and  r  by  x,  x',  and  r, 
respectively.  The  normal  derivative  of  the  Green’s  function  D  is  (Ref.  1) 


3D(r) 

da 


p-ikr  1 

— —  (ik  +  — )  cos  0, 
4ot  r 


(7) 


where  0  is  the  angle  between  the  normal  n  and  the  vector  r,  as  shown  in  Fig.  1. 

All  three  integral  equations  in  Eq.  4  are  needed  for  exterior  fluids.  The  surface 
equation  provides  the  fluid  impedance  at  the  fluid-structure  interface.  Since  the  surface 
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equation  exhibits  non-uniqueness  at  certain  discrete  characteristic  frequencies  (Ref.  9),  the 
interior  equation  is  used  to  provide  additional  constraint  equations  which  ensure  the  required 
uniqueness.  The  exterior  equation  is  used  to  compute  the  exterior  pressure  field  once  the 
surface  solution  (which  includes  the  fluid  pressure  and  its  gradient)  is  known. 

The  substitution  of  Eqs.  6  and  7  into  the  surface  equation  (4)  yields 

-  Js  p(x)  (ik  +  j)  cos  p  dS  =  i up  fs  vn(x)  -^-dS  +  pI?  x'  on  S.  (8) 

This  integral  equation  relates  the  total  pressure  p  and  normal  velocity  vn  on  S.  If  the 
integrals  in  Eq.  8  are  discretized  for  numerical  computation  (Ref.  1),  we  obtain  the  matrix 
equation  (for  the  exterior  fluid) 

Ep  =  Cvn  +pI;  (9) 

where  vector  p  (of  dimension  f  x  r)  is  the  vector  of  complex  amplitudes  of  the  total  pressure 
on  the  structure’s  wet  surface,  matrices  E  and  C  (both  f  x  f)  are  fully-populated,  complex, 
nonsymmetric,  and  frequency-dependent,  and  vector  pi  (f  x  r)  is  the  complex  amplitude  of 
the  incident  pressure  vector.  The  number  of  unknowns  in  this  system  is  f,  the  number  of 
wet  points  on  the  fluid-structure  interface. 

The  normal  velocities  v„  in  Eq.  9  are  related  to  the  total  velocities  v  by  the  same 
rectangular  transformation  matrix  G: 

vn  =  GTv,  (10) 

where  T  denotes  the  matrix  transpose.  If  velocities  v  and  v„  are  eliminated  from  Eqs.  1,  9, 
and  10,  the  resulting  equation  for  the  coupled  fluid-structure  system  is 

(E  +  CGtZ“1GA)  p  =  CG^F  +  Pl.  (11) 

This  equation  is  solved  for  the  total  surface  pressures  p,  since  the  rest  of  the  equation 
depends  only  on  the  geometry,  the  material  properties,  and  the  frequency.  Since  the  two 
right-hand  side  terms  in  Eq.  11  correspond  to  mechanical  and  incident  loadings,  only  one  of 
the  two  terms  would  ordinarily  be  present  for  a  given  case.  The  details  of  the  incident 
pressure  vector  pj  for  scattering  problems  were  presented  previously  (Ref.  2).  and  will  not  be 
repeated  here. 

The  velocity  vector  v  for  all  structural  DOF  is  recovered  by  solving  Eq.  1  for  v: 

v  =  Z^F-Z^GAp.  (12) 

The  surface  normal  velocity  vector  v„  is  recovered  by  substituting  this  solution  for  v  into  Eq. 

10. 


Modeling  Interior  Fluid 

The  theoretical  development  presented  in  the  preceding  section  can  be  modified  slightly 
to  account  also  for  an  interior  fluid.  The  wave  equation,  Eq.  3,  applies  also  to  interior  fluids. 
Although  all  three  integral  equations  in  Eq.  4  are  generally  needed  for  exterior  fluids,  only 
the  surface  equation  is  needed  to  represent  the  surface  impedance  of  interior  fluids.  Eq.  4a 
also  applies  to  interior  fluids  if  the  incident  pressure  pj  is  set  to  zero,  and  the  normal  vector 
n  is  negated.  That  is,  the  surface  integral  equation  applies  to  both  exterior  and  interior  fluids 
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so  long  as  the  unit  normal  is  always  directed  from  the  structure  into  the  fluid.  One  other 
consideration,  perhaps  unique  to  NASHUA,  is  that  wet  surface  curvatures  (which  are 
needed  in  the  calculation  of  the  "self"  terms  in  matrix  E)  are  negative  at  some  interior  points 
(Ref.  1). 

A  matrix  equation  similar  to  Eq.  (9)  is  therefore  obtained  for  the  interior  fluid  except 
that  the  incident  pressure  pi  is  zero.  The  fluid  matrices  E  and  C  are  different  for  exterior 
and  interior  domains  (even  if  the  separating  surface  S  has  infinitesimal  thickness)  because 
the  normals  are  of  opposite  sign. 


Two-Fluid  Formulation 


Denote  the  exterior  fluid  as  Fluid  1  and  the  interior  fluid  as  Fluid  2,  and  use  the 
subscripts  1  and  2  to  refer  to  these  two  domains.  Also  define  new  pressure  and  normal 
velocity  unknowns  p  and  vn  which  include  the  solutions  for  both  fluid  domains: 


Since  there  is  no  direct  fluid  coupling  between  the  interior  and  exterior  fluids,  and  the 
incident  pressure  vanishes  in  the  interior  domain,  Eqs.  1,  9,  10,  and  11  apply  also  to  the 
two-fluid  problem  if  the  new  definitions  in  Eq.  13  are  used,  and  the  matrices  A,  G,  E,  C, 
and  pi  are  re-defined  as 


A  = 


C  = 


Cx 


C\ 


(14) 


The  principle  benefit  of  formulating  the  two-fluid  problem  in  this  way  is  that  the  required 
modifications  to  extend  the  procedure  to  three  or  more  independent  fluid  domains  is  now 
clear. 


The  Far-Field  Calculation 

With  the  solution  for  the  total  pressures  and  velocities  on  the  surface,  the  exterior 
Helmholtz  integral  equation,  Eq.  4b,  can  be  integrated  to  obtain  the  radiated  (or  scattered) 
pressure  at  any  desired  location  x'  in  the  exterior  field.  We  first  substitute  Eqs.  5-7  into  Eq. 
4b  to  obtain 

1  e-ikr 

p(x')  =  f  [i w/>vn(x)  +  (ik  +  —  )p(x)  cos  /?]  —  dS,  x'  in  E.  (15) 

o  r  '4  /i  i 

In  applications,  however,  the  field  pressures  generally  of  interest  are  in  the  far-field,  so  we 
use  instead  the  asymptotic  form  of  Eq.  15  (Ref.  1): 
ikx' 

p(x')  =  — - — —  /  [/?cv„(x)  +  p(x)  cos  /?]e,kx  cos  0  dS,  x' in  E,  x' »  d,  (16) 

Attx  s 

where  d  is  a  characteristic  dimension  of  the  structure,  and  a  is  the  angle  between  the  vectors 
x  and  x'  (Fig.  1).  For  far-field  points,  cos  j3  is  computed  using  the  asymptotic  approximation 
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(17) 


a  X 

cos  B  — »■  n — -  . 

x 

For  both  Eqs.  15  and  16,  numerical  quadrature  is  used. 


OVERVIEW  OF  SOLUTION  PROCEDURE 

The  NASHUA  solution  procedure  uses  NASTRAN  to  generate  the  matrices  K,  M,  B, 
and  F  and  to  generate  sufficient  geometry  information  so  that  the  matrices  E,  C,  G,  A,  and 
Px  can  be  computed  by  a  separate  program  called  SURF.  Then,  NASTRAN  DMAP  is  used 
to  form  the  matrices  appearing  in  Eq.  11,  which  is  solved  for  the  total  pressures  p  (in  both 
fluid  domains)  using  the  block  solver  OCSOLV  (Ref.  10).  Next,  NASTRAN  DMAP  is  used 
to  recover  the  surface  normal  velocities  vn  and  the  vector  v  of  velocities  at  all  structural 
DOF  (NASTRAN’s  "g-set").  This  step  completes  the  surface  solution.  Then,  with  the  total 
pressures  and  velocities  on  the  (exterior)  surface,  the  asymptotic  (far-field)  form  of  the 
Helmholtz  exterior  integral  equation  is  integrated  in  program  FAROUT  to  compute  the  far- 
field  radiated  pressures.  Various  tables  and  graphical  displays  are  generated. 

The  overall  setup  of  the  solution  procedure  is  organized  into  four  steps.  In  Step  1,  a 
separate  NASTRAN  structural  model  is  prepared  and  run  for  each  unique  set  of  symmetry 
constraints  and  each  fluid  region.  Since,  for  general  three-dimensional  analysis,  up  to  three 
planes  of  reflective  symmetry  are  allowed,  there  would  be  one,  two,  four,  or  eight  such  runs 
for  each  fluid  region.  Since  the  purpose  of  this  step  is  to  generate  a  file  containing  geometry 
information  and  a  checkpoint  file  for  subsequent  use  in  the  other  steps,  the  only  difference 
between  the  two  runs  associated  with  a  given  symmetry  case  is  the  specification  of  the 
outwardly  directed  unit  pressure  load  which  defines  the  wet  surface  for  a  given  fluid  region. 

For  each  symmetry  case  and  drive  frequency,  several  programs  are  run  sequentially  to 
form  Step  2.  For  each  fluid  region,  the  SURF  program  reads  the  geometry  file  generated  by 
NASTRAN  in  Step  1  and,  using  the  Helmholtz  surface  and  interior  integral  equations, 
generates  the  fluid  matrices  Ej,  E 2,  Cj,  and  Co,  the  area  matrices  A]  and  A2,  the  structure- 
fluid  transformation  matrices  Gx  and  Go,  the  incident  pressure  vector  pn,  and  a  geometry 
file  to  be  used  later  by  the  far-field  integration  program  FAROUT  in  Step  3.  In  addition,  a 
partitioning  vector  is  generated  to  facilitate  the  merging  and  partitioning  of  the  various 
matrices  associated  with  the  two  fluid  domains. 

The  two  SURF  jobs  in  Step  2  are  followed  by  a  NASTRAN  job  which  takes  the 
structural  matrices  K,  M,  B,  and  F  from  Step  1  and  the  matrices  generated  by  the  SURF 
jobs  and  forms  the  matrices  in  Eq.  11,  where  the  definitions  in  Eq.  14  apply.  Eq.  11  is  then 
solved  for  the  total  surface  pressure  vector  p  by  program  OCSOLV,  which  is  a  general  out- 
of-core  block  solver  designed  specifically  for  large,  full,  complex,  nonsymmetric  systems  of 
linear,  algebraic  equations.  NASTRAN  is  then  re-entered  in  Step  2  with  p  so  that  the 
velocities  v  and  vn  can  be  recovered  using  DMAP  operations.  The  surface  pressures, 
normal  velocities,  and  full  g-set  displacements  are  then  reformatted,  sorted,  and  merged  into 
a  single  file  (for  each  symmetry  case)  using  program  MERGE.  Recall  that  there  are  one, 
two,  four,  or  eight  possible  symmetry  cases. 

Steps  1  and  2  are  repeated  for  each  symmetry  case.  After  all  symmetry  cases  have 
been  completed  and  merged,  program  FAROUT  (Step  3)  combines  the  symmetry  cases  and 
integrates  over  the  surface.  The  far-field  pressure  solution  is  obtained  by  integrating  the 
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surface  pressures  and  velocities  using  the  asymptotic  (far-field)  form  of  the  exterior 
Helmholtz  integral  equation,  Eq.  16.  Output  from  FAROUT  consists  of  both  tables  and  files 
suitable  for  various  types  of  plotting. 

The  remaining  steps  in  the  NASHUA  procedure  are  for  graphical  display.  Deformed 
structural  plots  of  the  frequency  response  are  obtained  by  restarting  NASTRAN  (Step  4) 
with  the  checkpoint  file  from  Step  1  and  a  results  file  from  FAROUT.  In  addition,  animated 
plots  can  be  generated  on  the  Evans  &  Sutherland  PS-330  graphics  terminal  using  the 
CANDI  program  written  for  the  DEC/VAX  computer  by  R.R.  Lipman  of  DTRC  (Ref.  11). 
X-Y  plots  of  various  quantities  (both  surface  and  far-field)  versus  frequency  may  be  obtained 
using  IPLOT  or  other  interactive  plotting  programs  (Ref.  12).  Polar  plots  of  the  far-field 
sound  pressure  levels  in  each  of  the  three  principal  coordinate  planes  can  also  be  generated 
using  the  interactive  graphics  program  FAFPLOT  (Ref.  1). 


DMAP  ALTERS 

Several  DMAP  alters  are  used  in  the  overall  NASHUA  procedure  to  implement  the 
precedure  described  in  preceding  section.  For  Step  1,  the  following  alter  is  used: 


ALTER 

ALTER 

ALTER 

GP3 

ALTER 

SSG1 

SSG2 

OUTPUT2 

OUTPUT2 

OUTPUT2 

PARAMR 

COND 

PARAMR 

DIAGONAL 

ADD 

RBMG2 

SSG3 

SDR1 

TA1 

DSMGl 

EQUIV 

COND 

MCE2 

LABEL 

EQUIV 

COND 


1  S  NASHUA  STEP  1,  COSMIC  1988  RF8  (REVISED  12/7/89) 

2,2  S  DELETE  PRECHK 
21,21  $  REPLACE  GPS 

GEOM3,EQEXIN,GEOM2/SLT,GPTT/S,N,NOGRAV/NEVER=l  $ 
117,117  S  REPLACE  FRRD 

SLT,BGPDT,CSTM,SIL,EST,MPT,GPTT,EDT,MGG,CASECC,DIT,/ 
PG„„/LUSET/NSKIP  S  PG 
USET,GM,YS,KFS,GO,DM,PG/QR,PO,PS,PL  S  PL 
AXIC,BGPDT,EQEXIN,USET,PG  S 
PL,CSTM,ECT„  $ 

$ 

//*EQV/C,Y,IISP=0./0.////NOHSP  $ 

LBL4D,NOIISP  $  SKIP  DIFF.  STIFF.  IF  NO  HYDRO.  P 
//*COMPLEX//C,Y,HSP=0./0./IiSPC  S  IISP+PO 
KAA/KDIAG/*SQUARE*/1.0  $ 
KAA,KDIAG/KAAD/(1.0,0.0)/(l.E-6,0.)  $ 

KAAD/LLL  $  FACTOR  KAA 

LLL,KAAD,PL,LOO,KOO,PO/ULV,UOOV,RULV,RUOV/OMlT/ 
V,Y,IRES=-1/1/S,N,EPSI  $  STATIC  SOLUTION 
USET,PG,ULV,UOOV,YS,GO,GM,PS,KFS,KSS,/UGV,PGG,QG/l/ 
*BKL0*  S  RECOVER  DEPENDENT  DISPLACEMENTS 
ECT,EPT,BGPDT,SIL,GP'1T,CSTM/X1,X2,X3,ECPT,GPCT/LUSET/ 
NOSIMP/O/NOGENL/GENEL  S  TABLES  FOR  DIFF.  STIFF. 
CASECC,GPTT,SIL,EDT,UGV,CSTM,MPT,ECPT,GPCT,D1T/KDGG/ 
S,N,DSCOSET  $  DIFF.  STIFF.  MATRIX 

KDGG.KDNN/MPCF2  /  MGG,MNN/MPCF2  S  EQUIV  IF  NO  MPC’S 
LBL1D,MPCF2  S  TRANSFER  IF  NO  MPC’S 
USET,GM,KDGG,„/KDNN,„  $  MPC’S  ON  DIFF.  STIFF. 

LBL1D  S 

K DN N , KDFF/ SI NGLE  /  M NN , M FF/SIN GLE  $  EQUIV.  IF  NO  SPC’S 
LBL2D, SINGLE  $  TRANSFER  IF  NO  SPC’S 
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SCE1  USET,KDNN,„/KDFF,KDFS,KDSS,„  $  SPC’S  AND  DIFF.  STIFF. 

LABEL  LBL2D  $ 

EQUIV  KDFF.KDAA/OMIT  /  MFF,MAA/OMIT  $  EQUIV.  IF  NO  OMITS 

COND  LBL3D,OMIT  $  TRANSFER  IF  NO  OMITS 

SMP2  USET,GO,KDFF/KDAA  $  OMITS  AND  DIFF.  STIFF. 

LABEL  LBL3D  $ 

PARAMR  //*SUBC*////MHSPC//HSPC  $  NEGATE  HYDRO.  P 

ADD  KDD,KDAA/NEWKDD/(1.0,0.0)/MHSPC  S 

ADD  KFS,KDFS/NEWKFS/(1.0,0.0)/MHSPC  $ 

EQUIV  NEWKDD,KDD//NEWKFS,KFS  $ 

LABEL  LBL4D  $  END  OF  DIFF.  STIFF.  EFFECTS  (HSP) 

DIAGONAL  KDD/IDENT/*SQUARE*/0.  $  D-SET  IDENTITY 

ADD  IDENT,/IDM/(1.0,0.0)  $  ANOTHER  D-SET  IDENTITY 

ADD  IDENT,/ZERO/ (0.0, 0.0)  $  D-SET  ZERO  MATRIX 

FRRD  CASEXX,USETD,DLT,FRL,GMD,GOD,IDENT,ZERO,IDM„DIT/ 

UDVF,PSF,PDF,PPF/*DISP*/*DIRECT*/LUSETD/MPCF1/ 
SINGLE/OMIT/NONCUP/FRQSET  $  PDF,  KDD=I,  BDD=0,  MDD=i 
CHKPNT  MDD,KDD,BDD,PDF,PSF,PPF,EQDYN,USETD,GOD,GMD  $ 

CHKPNT  KFS,BGPDT,ECT,EQEXIN,GPECT,SIL  S 

EXITS 

END  ALTER  $ 


The  above  alter  does  not  depend  on  whether  the  fluid  is  interior  or  exterior  to  the  structure. 
The  Step  2  alters,  however,  depend  on  whether  an  interior  fluid  is  present.  For  Step  2A,  the 
following  alter  is  used: 

ALTER  1  S  NASHUA  STEP  2A,  COSMIC  1988  RF8  (REVISED  11/7/89) 

ALTER  2,167  S  REPLACE  ALL  AFTER  ’BEGIN’  AND  BEFORE  ’END’ 
INPUTT2  /DAT2„„//13  S  INTERNAL  FLUID 

INPUTT2  /DAT„„//11  $  READ  SURF  MATRIX  FROM  UT1 

MATPRN  DAT,DAT2,„  $ 

PARAML  DAT//*DMP/1/S/RIGD  $  GET  RIGID  FLAG 

PARAMR  / /*EQV /RIGD/0. /// /ELAST  $  SET  ELAST=-1  IF  ELASTIC 

COND  LBL9D, ELAST  $  IF  ELASTIC,  JUMP  OVER  RIGID/SOFT 

PARAMR  //*EQ*//RIGD/2.////SOFT  $  SET  SOFP=-l  IF  SOFT  BD. 

COND  LBL9E.SOFT  $  IF  SOFT  BOUNDARY,  JUMP  OVER  RIGID 

INPUTT2  /E,PI,VEKC,,//11  S  READ  SURF  MATRICES  FROM  UT1 
OUTPUT2  PI,E,„  H-1  $  INPUTT2  FILE  TS  OVER-WRI1TEN  (UT1) 

OUTPUT2  „„  //-9  S  EOF 

CHKPNT  DAT,VEKC  S 

EXIT  $ 

LABEL  LBL9E  $  BEGINNING  OF  SOFT  ANALYSIS 

INPUTT2  /CT,PI,VEKC„//11  S  READ  SURF  MATRICES  FROM  UT1 

TRNSP  CT/C  S 

ADD  PI,/MPI/ (-1 .0,0.0)  S  NEGATE  PI 

OUTPUT2  MPI,C„,  //-l  S  INPUTS  FILE  IS  OVER-WRITTEN  (UT1) 

OUTPUT2  „„  //-9  $  EOF 

CHKPNT  DAT.VEKC  $ 

EXIT  S 

LABEL  LBL9D  $  BEGINNING  OF  ELASTIC  ANALYSIS 
INPUTT5  /G2,A2,„//14  S  INTERNAL  FLUID 
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INPUTT2  /C2T,E2,PI2,V£KC2,//13  $  INTERNAL  FLUID 

INPUTT5  /Gl,Al„,//12  $  READ  SURF  MATRICES  FROM  UT2 

INPUTT2  /ClT,El,PIl,VEKC,FVEC//ll  $  READ  SURF  MATRICES 

MATPRN  FVEC,„,  $ 

MERGE  A1,„A2,FVEC,/A/-1  $ 

MERGE  E1,„E2,FVEC,/E/-1  $ 

MERGE  C1T,,,C2T,FVEC,/CT/-1  S 

MERGE  Gl,,G2„FVEC,/G/0  S 

MERGE  PIl„„,FVEC/PI/0  $ 

MATPRN  VECM,VECS,„  S 

PARAML  DAT//*DMP/l/2/FREQ  S  GET  FREQ  FROM  DAT 

PARAMR  //*COMPLEX*//FREQ/0./FREQC  $  FREQ+PO 

PARAMR  //*MPY C*// // W/FREQ C/ (6 .283185,0.)  $  OMEGA 

PARAMR  //*MPYC*////IW/W/(0.,1.)  $  POMEGA 

PARAMR  //*MPYC!7///WW/W/W  $  OMEGA**2 

PARAMR  //* SUBC*// //MWW//WW  S  -OMEGA**2 

A.DD5  MDD,KDD,BDD„/Y/MWW/(1.0,0.0)/IW  $ 

MPYAD  G,A,/GA/0  $ 

DECOMP  Y/L,U/1//S,N,MINDIAG///S,N,SING  $ 

FBS  L,U,GA/YTGA/1  $ 

FBS  L,U,PDF/YIF/1  $ 

ADD  YIGA,/ZIGA/IW  S 

ADD  YIF,/ZIF/IW  $ 

MPYAD  G,ZIGA,/GTZIGA/1  $ 

MPYAD  CT,GTZIGA,E/II/1  $  LIIS 

MPYAD  G,ZIF,/GTZIF/1  S 

MPYAD  CT,GTZIF,/Q/1  $  MECHANICAL  RIIS 

MERGE  DUM„PDF„VECM,/PDF1/1  $  MERGE  IN  0  COLUMNS 

MERGE  DUM„PSF„VECM,/PSF1/1  S  MERGE  IN  0  COLUMNS 

MERGE  DUM„PPF„VECM,/PPF1/1  $  MERGE  IN  0  COLUMNS 

EQUIV  PDF1,PDF//PSF1,PSF//PPF1,PPF  $ 

MERGE  DUM„Q,,VECM,/RIIS1/1  S  MERGE  IN  ZERO  COLUMNS 

MERGE  DUM„GTZIF„VECM,/GTZIFE/1  $  MERGE  IN  0  COLUMNS 
MERGE  DUM„PI„VECS,/RHS2/1  S  MERGE  IN  ZERO  COLUMNS 

ADD  RHS1,RHS2/RHS  S  ADD  MECII.  AND  INC.  RHS 

EQUIV  USETD,DUMI//GOD,DUM2//GMD,DUM3//KFS,DUM4  $ 

OUTPUT2  RHS, II,,,  //-l  $  INPU'IT2  FILE  IS  OVER-WRI1TEN  (UTI) 

OUTPUT2  „„  //-9  S  EOF 

CIIKPNT  GTZIGA,GTZIFE,GA,PDF,L,U,PSF,DAT,VEKC,FVEC  S 

CHKPNT  USETD,GOD,GMD,KFS  $ 

ENDALTER  $ 

The  differences  between  this  alter  and  one  used  for  submerged  evacuated  structures  are  due 
to  the  need  to  read  and  combine  two  sets  of  SURF  matrices,  one  for  each  fluid  domain.  For 
Step  2B,  the  following  alter  is  used: 

ALTER  1  S  NASHUA  STEP  2B,  COSMIC  19SS  RF8  (REVISED  11/7/89) 

ALTER  2,167  $  REPLACE  ALL  AFrER  'BEGIN’  AND  BEFORE  'END' 

INPUTT2  /PC„„//11  $  READ  PRESSURES  FROM  BLOCK  SOLVER  (UTI) 

PARTN  PC„FVEC/Pl„,/0  $  REMOVE  INTERNAL  FLUID  DOF 
PARTN  Pl„VEKC/P„,/0  $  REMOVE  CHIEF  DOF  FROM  P 
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COND  LBL9D,ELAST  S  IF  ELASTIC,  JUMP  OVER  RIGID/SOFT 

OUTPUT2  DAT,P,„  ll-l  $  INPUTT2  FILE  IS  OVER-WRITTEN  (UT1) 

OUTPUT2  „„  //- 9  $  EOF 

MATPRN  DAT,P,„  S  FOR  SOFT  BOUNDARY,  P  REPRESENTS  VELOCITY 
EXIT  $ 

LABEL  LBL9D  $  ELASTIC  ANALYSIS 

MPYAD  GTZIGA,PC,GTZIFE/VNC/0/-l  S  NORMAL  VELOCITIES 

MPYAD  GA,PC,PDF/FA/0/-l  $  A-SET  FORCES 

FBS  L,U,FA/UDVF/1  $  A-SET  DISPLACEMENTS 

SDR1  USETD„UDVF,„GOD,GMD,PSF,KFS„/UPVC„QPC/l/ 

^DYNAMICS*  $ 

PARTN  VNC,,FVEC/Vl„,/0  S  REMOVE  INTERNAL  FLUID  DOF 

PARTN  Vl„VEKC/VN„,/0  S  REMOVE  CHIEF  DOF  FROM  VN 

OUTPUT2  DAT,P,VN,UPVC,  ll-l  $  INPUTT2  FILE  IS  OVER- WRITTEN 

OUTPUT2  „„  //- 9  $  EOF 

MATPRN  DAT,P,VN„  $ 

ENDALTER  S 

This  alter  differs  from  one  for  evacuated  structures  because  of  the  presence  of  several  matrix 
partitionings  to  remove  the  internal  fluid  DOF  from  the  solution  vectors  before  the  solutions 
are  merged  with  the  results  for  other  frequencies. 


NUMERICAL  EXAMPLE 

Here  we  illustrate  and  validate  the  two-fluid  boundary  element  formulation  developed 
above  by  solving  the  problem  of  acoustic  scattering  from  a  submerged  fluid-filled  spherical 
thin  shell.  The  incident  loading  is  a  time-harmonic  planar  wavclrain,  as  shown  in  Fig.  2. 
The  specific  problem  solved  has  the  following  characteristics: 


shell  mean  radius  (a) 

5  m 

shell  thickness  (h) 

0.15  m 

shell  Young’s  modulus  (E) 

2.07  x  10u  N/nr 

shell  Poisson’s  ratio  (//) 

0.3 

shell  density  ( ps ) 

7669  kg/m3 

shell  loss  factor  (tj) 

0.01 

fluid  density  (p) 

1000  kg/m3 

fluid  sound  speed  (c) 

1524  m/s 

The  same  fluid  is  used  for  both  the  exterior  and  interior  fluid  domains.  The  solution  of  this 
problem  exhibits  rotational  symmetry  about  the  spherical  axis  parallel  to  the  direction  of 
wave  propagation.  The  benchmaik  solution  to  which  the  numerical  results  will  be  compared 
is  a  series  solution,  the  derivation  of  which  is  summarized  in  the  next  section. 


Series  Solution 

The  series  solution  for  scattering  from  a  submerged  evacuated  spherical  thin  shell  was 
presented  by  Junger  and  Feit  (Ref.  13).  Here  we  summarize  that  solution  and  indicate  the 
modification  necessary  to  include  the  addition  of  an  interior  fluid  which  fills  the  spherical 
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Fig.  2.  Plane  Wave  Scattering  from  a  Fluid-Filled  Spherical  Shell 

volume. 


In  general,  the  series  solution  for  plane  wave  scattering  from  a  submerged,  evacuated, 
spherical  thin  shell  involves  computing  the  impedances  of  the  shell  and  exterior  fluid,  the 
scattered  field  due  to  rigid  body  effects,  and  the  radiated  field  due  to  elastic  (vibrational) 
effects.  The  shell  impedance  (the  ratio  of  pressure  to  normal  velocity)  for  the  nth 
axisymmetric  shell  mode  is 


ifec„  h  [tf-fflW)2]  [tf-(n!?>)2] 
a  [O'— (l+/?~)(i'-}-Xn— 1] 


where  ps  is  the  structural  mass  density,  cp  =vE/[ft(l-/^)]  ,  E  is  Young’s  modulus,  v  is 
Poisson’s  ratio,  Q=  wa/cp  is  th£_ dimensionless  frequency,  h  is  the  shell  thickness,  a  is  the 
shell  mean  radius,  /?  =  h/(a  m.2),  and  X„=n(n+1).  The  quantities  and  ftp  are  the 
upper  and  lower  shell  resonance  dimensionless  frequencies,  respectively.  They  are  the 
solutions  of  the  characteristic  equation 


17'  —  [l+3t/+X„— /?2(1—  y— X„—  Ai]  fi2 

+  (Xn  — 2)  ( 1  —ir>  [Xn  3  — 4Xn  2  H-X(5— z/2 )  — 2  ( 1  —\r )]  =  0.  (19) 


The  impedance  of  the  exterior  fluid,  found  by  using  the  Green’s  function  and  identity  for  the 
exterior  fluid,  is 


z 


n 


i/7C 


w 

h'n(ka)  ’ 


(20) 


where  hn  is  the  Bessel’s  function  of  the  third  kind  of  order  n. 


Thus,  Junger  and  Feit  showed  that  the  far-field  scattered  pressure  is 
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pOM)= 


,  R>>a,  (21) 


ielkRp0  (2n+l)Pn(cos0) 

k&  n=0  h'n(ka) 


j'n(ka) 


_ Pc  ■ 

(Zn+Zn)(ka)2h'n(ka) 


where  R  is  the  distance  to  the  field  point,  0  is  the  angle  from  the  z-axis,  p0  is  the  incident 
pressure,  Pn  is  the  Legendre  polynomial  of  order  n,  and  jn  is  the  Bessel’s  function  of  the  first 
kind  of  order  n.  The  two  terms  in  the  bracketed  expression  correspond  to  rigid  body  and 
radiated  effects,  respectively. 


The  above  expression  for  the  pressure  scattered  from  an  evacuated  shell  can  be 
extended  to  include  the  effects  of  the  interior  fluid  merely  by  replacing  the  exterior  fluid 
impedance  zn  in  Eq.  21  with  the  sum  of  the  fluid  impedances  for  the  exterior  and  interior 
fluids.  It  can  be  shown,  by  using  the  Green’s  function  and  identity  for  the  interior  domain, 
that  the  interior  impedance,  denoted  fn,  is  given  by 


?n  =  pc 


jnCka) 
j'n(ka)  ’ 


(22) 


We  ;,ote  the  resemblance  between  Eqs.  20  and  22  for  the  exterior  and  interior  domains, 
respectively. 


The  computer  program  used  to  evaluate  this  series  solution  is  a  modification  of  a 
program  called  SCATSPIIERE  written  by  F.M.  Henderson,  a  retired  employee  of  DTRC. 
SCATSPIIERE  in  turn  is  a  variant  of  an  earlier  program  called  RADSPIIERE  (Ref.  14)  for 
computing  the  radiation  from  an  internally-driven  submerged  spherical  shell. 


Numerical  Solution 

A  NASTRAN  finite  clement  model  of  the  spherical  shell  was  prepared  using  40 
axisymmetric  conical  shell  elements  spanning  the  180  degrees  between  the  two  poles  of  the 
sphere.  Due  to  the  axisymmetry  of  the  incident  pressure  loading,  only  the  N  ==  0  harmonic 
was  required.  Since  all  structural  points  are  in  contact  with  both  interior  and  exterior  fluids, 
the  resulting  model  therefore  had  205  independent  structural  degrees  of  freedom  (DOF)  and 
41  fluid  DOF  for  each  of  the  two  fluid  domains.  System  matrices  for  the  exterior  fluid  were 
also  augmented  by  the  addition  of  four  constraint  equations  associated  with  interior  Chief 
points  to  ensure  uniqueness  of  the  integral  representation  at  the  upper  frequencies.  The 
nondimcnsional  frequency  range  0<ka<5  was  swept  using  a  frequency  increment  of  about  ka 
=  0.05  with  NASHUA  and  ka  =  0.005  with  the  scries  solution.  Since  the  series  solution  is 
converged,  we  treat  it  as  an  "exact"  solution  for  this  problem. 

The  comparison  between  the  computed  and  exact  solutions  is  presented  is  Figs.  3  and 
4,  which  plot  the  frequency  response  of  the  nondimensional  scattered  pressure  pr/(pQa), 
where  p  is  the  far-ficld  scattered  pressure  at  distance  r  from  the  origin,  p0  is  the  incident 
pressure,  and  a  is  the  mean  radius  of  the  spherical  shell.  These  two  figures  show  very  good 
agreement  between  the  two  scattering  solutions  in  the  backward  ( 0  =  0)  and  forward  (0  =  180 
degrees)  directions.  In  fact,  the  computed  and  series  solutions  are  virtually  indistinguishable 
from  each  other. 
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Fig.  3.  Forward  Scattering  from  a  Fluid-Filled  Spherical  Shell 


Fig.  4.  Backward  Scattering  from  a  Fluid-Filled  Spherical  Shell 

DISCUSSION 

A  very  general  computational  capability  has  been  described  for  predicting  the  sound 
pressure  field  radiated  or  scattered  by  arbitrary,  submerged,  fluid-filled,  three-dimensional 
elastic  structures  subjected  to  time-harmonic  loads.  The  structure  is  modeled  with 
NASTRAN  (in  all  the  generality  that  NASTRAN  allows)  in  combination  with  boundary 
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element  models  of  both  interior  and  exterior  fluid  domains.  Sufficient  automation  is 
provided  so  that,  for  many  structures  of  practical  interest,  an  existing  structural  model  can 
be  adapted  for  NASHUA  acoustic  analysis  within  a  few  hours. 

One  of  the  many  benefits  of  having  NASHUA  linked  with  NASTRAN  is  the  ability  to 
integrate  the  acoustic  analysis  of  a  structure  with  other  dynamic  analyses.  Thus  the  same 
finite  element  model  can  be  used  for  modal  analysis,  frequency  response  analysis,  linear 
shock  analysis,  and  underwater  acoustic  analysis.  In  addition,  many  of  the  pre-  and 
postprocessors  developed  for  use  with  NASTRAN  become  available  for  NASHUA  as  well. 
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MONITORING  OF  RITZ  MODAL  GENERATION 


Mladen  Chargin  and  Thomas  G.  Butler 

NASA  Ames  Research  Center  Butler  Analyses 


In  applying  a  Ritz  modal  expansion  to  the  solution  of  a 
transient  response,  there  is  a  problem  as  to  how  many  modes  are 
needed  to  obtain  accuracy  to  within  a  specified  percentage. 

One  of  us,  Chargin,  has  suggested  a  method  based  on  the 
characteristics  of  the  forcing  function.  The  method  can  be 
incorporated  into  the  Ritz  generation  algorithm  such  that  it  will 
automatically  monitor,  regulate  and  terminate  the  process 
according  to  a  specified  tolerance. 

FORCING  FUNCTION  CHARACTERISTICS 

Assume  that  the  forcing  function  F(x,t)  can  be  repre¬ 
sented  as  a  product  of  a  spatial  function  and  a  temporal  func¬ 
tion;  i.e. 

(1)  F(x,t)  =  F(x)  .  f (t) . 

Develop  a  criterion  based  upon  measuring  the  amount  of  power 
developed  in  the  forcing  function  F(x,t).  The  total  power  in 
F(x,t)  is  the  product  of  the  "power"  in  F(x)  and  f(t)  owing  to 
the  assumption  in  equation  (1).  The  scheme  in  outline  is  to 
measure  the  temporal  power  and  the  spatial  power  in  F(x,t)  separ¬ 
ately  then  compare  the  corresponding  power  in  the  Ritz  modes 
against  power  in  each  component  of  the  forcing.  Generation  of 
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additional  Ritz  modes  will  continue  until  the  power  criteria  are 
met.  First  a  measure  is  taken  of  the  total  power  in  the  temporal 
function. 


Temporal  Power 


(2) 


P(T)  = 


[m,] 


dt 


where  t  is  the  interval  over  which  the  transient  will  act.  There 
could  very  well  be  a  separate  temporal  function  for  each  spatial 
function.  In  that  case  of  multiple  loadings  P(t)  of  equation  (2) 
would  be  a  vector  "1"  long. 


In  order  to  tailor  this  power  to  our  use  as  a  guide  in 
selecting  Ritz  modes  it  will  be  useful  to  measure  the  amount  of 
temporal  power  as  a  function  of  frequency.  Expand  the  temporal 
function  in  a  Fourier  Series  and  sum  the  power  versus  the  expan¬ 
sion  multiple: 


(3) 


f(t)  =  aQ  +  a-^cos  wt  +  b-^sin  wt  +  a2Cos  2wt  +  b2Sin  2wt 


+ . +a  cos  wt  +  b  sin  wt  . 

n  n 


The  power  within  a  band  0  to  nw  is: 

n  ,  x  2 

(4) 


P  =  >  fa. cos  iwt  +  b.sin  iwt)  =  >  fa?  +  b?)  . 
n  |=0l  i  i  )  I=0l  i  iJ 


One  can  compare  the  amount  of  power  within  a  given  band  P^  with 
the  total  power  P(t).  When  Pn  is  within  close  range  of  P(x),  say 
1%,  then  the  analyst  can  be  satisfied  that  the  frequency  range  of 
the  truncated  temporal  forcing  function  is  sufficiently  broad  to 
encompass  the  temporal  requirements  of  the  forcing. 
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For  a  certain  frequency  f  =  5^  ;  kv2t  =  .99.  Let  this  fre- 

J  n  2ir  P(t) 

quency  be  r  . 


Spatial  Power 


Now  turn  to  the  spatial  distribution  of  the  forcing  function 
F(x),  and  develop  a  measure  that  is  called  spatial  power.  Make 
an  additional  assumption  that  all  points  being  loaded  have  mass. 


(5) 


CH(x)3  =  CF( x) 3  CM3  CF(x) 3 


where  |CjaJ  =  jVfxJ^Wj  is  a  coefficient  matrix  that  will  be  shown 

to  be  usefual  later.  H(x)  is  a  matrix  if  there  are  a  number  of 
loading  cases  "1". 
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At  this  point  we  have  2  measures  of  the  forcing  fuction. 
We  have  f^  reflecting  the  desired  frequency  content  and  IT(x) 
reflecting  the  desired  level  of  power  to  activate  the  structural 
mass.  Now  it  is  time  to  test  the  adequacy  of  the  number  of  Ritz 
modes  to  repond  to  the  forcing  at  these  power  levels. 


First,  get  a  spatial  measure  of  the  Ritz  modes  0^.  A 
simple  scheme  is  to  pattern  the  measure  after  equation  (5)  but 
substitute  the  matrix  of  k  spatial  Ritz  functions  for  the  post- 
multiply  operation  instead  of  the  spatial  component  of  loading 
F(x) . 


(6) 


CF(x) JTCM3 


k=l 


The  way  to  use  R^  is  to  compare  each  diagonal  term  of  II(x)  with 
the  corresponding  diagonal  term  of  Rr  to  zee  if  the  ratio  is 
within  a  specified  tolerance;  i.e. 


(7) 


<  e  for  every  1. 


Keep  generating  additional  Ritz  modes  0^  k  >  r  until  an  r  has 
been  reached  for  which  every  diagonal  term  satisfies  the 
criterion. 


Monitoring  Function 


If  the  0^  satisfy  the  spatial  power  criterion  of  F(x,t)  it  does 
not  necessarily  hold  that  0^  will  simultaneously  satisfy  the 
temporal  requirements  of  F(x,t).  Therefore,  use  the  0^  which 
have  met  the  spatial  power  requirements  and  obtain  an  estimate  of 
its  frequency  content  by  setting  up  the  frequency  equation.  (Jse 


168 


MONITORING  OF  RITZ  MODAL  GENERATION 


an  estimator  that  is  much  less  demanding  than  solving  che 
eigenvalue  problem.  Merely  do  the  following. 

t  ll 

Compute  a  k —  order  Ritz  generalized  mass  and  stiffness;  i.e. 


(81  M  ■  KMeJ  and  h] ' 


and  construct  a  test  matrix  called 


threshold  frequency,  f  ,  determined  from  the  temporal  power  P^, 


J 


which  involves  the 


-  (2irf  )\1  = 
o  kj 

Decompose  IS^  and  extract  the  value  of  the  output  parameter 
NBRCHG  issued  by  the  DECOMP  module.  Param  NBRCHG  reports  the 
number  of  negative  values  on  the  factor  diagonal  of  which  is 
tantamount  to  the  number  of  sign  changes  or  zero  crossings  in  the 
characteristic  equation.  If  the  value  of  NBRCHG  =  k,  it 
indicates  that  the  frequencies  of  all  k  Ritz  modes  are  less  than 
fQ.  One  would  be  inclined  to  want  the  frequency  content  of  the 
Ritz  modes  to  bracket  f  ;  i.e.  some  modes  with  a  frequecy  f  . 
This  implies  that  one  would  seek  to  have  the  value  of  NBRCHG  to 
be  less  than  the  order  of  matrix  S^.  This  idea  can  be  built  into 
the  Ritz  generation  routine  as  a  test  as  to  whether  enough  modes 
have  been  generated  to  within  a  certain  margin  a  such  that  f  > 
afQ,  where  the  user  specifies  a. 


one  would 


not  want  to  repeat  the  eigenvalue 

because 


ekK  ek 


Obviously, 

estimate  each  time  a  new  Ritz  mode  is  obtained, 

T 

and  9^M  6^  could  become  expensive  as  k  becomes  large.  Therefore 

develop  a  scheme  whereby  the  spatial  power  is  reduced  by  some 

factor  v  >  1;  i.e.  e  =  e  , ,/y.  Typically,  one  can  use  a  value 
1  new  old  '  J 
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of  2  to  4  for  y*  Once  the  spatial  error  e  ,  is  satisried, 

L  new 

repeat  the  eigenvalue  estimate  test. 

Conclusion 


A  scheme  has  been  proposed  to  monitor  che  adequacy  of  a 
set  of  Ritz  modes  to  represent  a  solution  by  comparing  the 
quantity  generated  with  certain  properties  involving  the  forcing 
function.  In  so  doing  an  attempt  has  been  made  to  keep  this 
algorithm  lean  and  efficient,  so  that  it  will  be  economical  to 
apply.  Using  this  monitoring  scheme  during  Ritz  Mode  generation 
will  automatically  ensure  that  the  k  Ritz  modes  9^  that  are 
generated  are  adequate  to  represent  both  the  spatial  and  temporal 
behavior  of  the  structure  when  forced  under  the  given  transient 
condition  defined  by  F(x,t). 
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